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I. INTRODUCTION 



A. BACKGROUND 

Extensive work has been accomplished recently in the area of modeling and 
predicting quarterly overhead costs for aircraft manufacturers. Overhead costs are 
generally predicted utilizing estimated overhead rates applied to labor hours or costs 
over functional categories such as manufacturing or engineering. Changes in operating 
rates cause overhead rate changes which may be observed only after a significant lag in 
time. Recent work done in this area has been to estimate total overhead as a function 
of the number of direct manufacturing personnel [Ref. 1]. 

Data collected for purposes of estimating overhead costs, since it is of a time 
series nature, can be expected to exhibit some degree of autocorrelation. Specifically, 
data collected for recent purposes has exhibited first order autoregressive AR(1), fourth 
order autoregressive AR(4), or a combination of these processes. An AR(1) process 
occurs when the errors in adjacent time periods are related. This type of relationship 
would be expected in yearly cost and operating data. A special form of the fourth 
order autoregressive would be expected in data of a seasonal nature, i.e. quarterly 
observations. Errors in this special case of the AR(4) process would be related to 
errors which occured four quarters previously. 

The method of ordinary least squares (OLS) is used to regress the data and 
estimate overhead costs. Ordinary least squares' validity is based on certain key 
assumptions: 

• Errors are distributed independently of the explanatory' variables with zero mean 
and constant variance. 

• Successive errors are independent of each other. 

With autocorrelation present, assumption two is violated and the ordinary least 
squares procedure breaks down at three points: 

• Estimates of regression coefficients, though unbiased are not efficient. 

• The usual formula for the variance of an estimate no longer applies and is liable 
to seriously underestimate the true variance. 

• t and F distributions, used for making confidence statements are no longer valid 
(Ref. 2). 

Since the OLS procedure breaks-down with the presence of autocorrelation, it is 
important that its presence be detected and the disturbance in the data subsequently 
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corrected. Without accounting for these disturbances in the data, estimates of the 
regression coefficients will not have minimum variance. 

The Durbin-Watson statistic is utilized to detect the presence of AR(1) behavior. 
Additionally, it has been generalized to detect higher order autoregressive processes. 
Durbin and Watson's work was based on the findings of T.W. Anderson which showed 
that the statistic 

e'Ae / e'e, (eqn 1.1) 

where e is the column vector of residuals from a regression and A is a certain real 
symmetric matrix, provides a test that is uniformly most powerful against certain 
alternative hypotheses [Ref. 2). 

The general model upon which Durbin and Watson's work was based is 



where X t is a (1 x K) non-stochastic matrix, P is a (K x 1) vector and i is one for an 
AR(1) process and four for an AR(4) process or a mixture of both processes. The 
error component, £ t , is specified by the form above and r| t is distributed with zero 
mean and constant variance. 

Typically an ordinary least squares regression is performed on selected 
independent variables. If the presence of autocorrelation is suspected, the data may be 
tested for a particular AR process utilizing the appropriate form of the Durbin-Watson 
statistic. For an AR(1) process the form of the Durbin-Watson statistic is 




(eqn 1.2) 



where 



£ t = 0 t c t _ i + n t . t= 1, T 



(eqn 1.3) 




(eqn 1.4) 



where 




(eqn 1.5) 



and 
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(eqn 1.6) 




The estimator of P from the ordinary least squares regression is $ The null hypothesis 
for the test is that of zero autocorrelation in the residuals against an alternative 
hypothesis that a first order autoregressive process exists. 

If the presence of an autoregressive process is detected in the regression residuals, 
the linear model must be reestimated after transforming the original data. For an 
AR(1) process, the independent variables are transformed as 



X," - X t (l - p, 2 ) 1 / 2 t-1 



(eqn 1.7) 



and 

X t = Xj - t= 2,...T. (eqn 1.8) 

In a similar fashion, the dependent variables are transformed as 

Y t * = Y t (l - Pl 2 ) 1/2 t=l (eqn 1.9) 

and 

Y t * = Y t - PjY t _! t = 2,...T. (eqn 1.10) 

In order to accomplish this transformation an estimate for pj is required. 
Though this parameter can be estimated in several ways, the most popular method due 
to simplicity of calculation is 

Pj = 1 - .5dj, (eqn 1. 1 1) 

where dj is given in equation 1.4 . In addition to the ease of calculation, this estimator 
performs as well as more complex estimators which are available |Ref. 3]. 

After the data is transformed, generalized least squares(GLS) is used to 
reestimate model parameters. After GLS is performed, the model is rechecked for the 
presence of autocorrelation using the Durbin-Watson statistic. 
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B. FOCUS AND SCOPE OF RESEARCH 

As previously mentioned, regression data is at times influenced by multiple 
autoregressive processes. In overhead cost estimating, both an AR(1) and AR(4) 
process have been observed to simultaneously influence certain data sets [Ref. 1], In 
the linear least squares model 



Y t = X t p + £ t , (eqn 1.12) 

suppose the errors are serially correlated and of the form 

£ t = 0 l £ t-l + 6 4 £ t-4 + n t > ( ec l n 113 ) 

where the effects of the second and third prior quarters are negligible compared to the 

effect of the most recent and year-earlier quarters. The rj t 's in the model, as in the 
previous version, are independent and distributed with zero mean and constant 
variance. 

It is advised [Ref. 4: p. 211] that if this special form of the AR(4) process exists, 
that it be tested for by a two-step procedure. In this procedure, after initially 
performing a OLS regression, the residuals are tested for an AR(1) disturbance utilizing 
the Durbin-Watson statistic for a first order process. If the influence of an AR(1) 
process is detected, a value for Pj is estimated and the original data is transformed. 
GLS is then performed on the transformed data and the residuals are tested for the 
influence of a fourth order autoregressive process by means of Wallis' test [Ref. 5]. If 
an AR(4) process is present, a value for p^ is estimated and the data is transformed a 
second time. At this point, the procedure is repeated to determine if either form of 
autocorrelation is still present. 

Unless the individual performing the test has a priori knowledge that both an 
AR(1) and AR(4) process exist within the data, it is very likely that one or the other 
would be overlooked. Additionally, peculiarities in the data may preclude the detection 
of one or the other processes. If this occurred, the GLS regression estimates as 
previously discussed would be inefFicient. In order to avoid the problem above, it 
would be desirable to have a procedure to detect simultaneously both an AR(1) and 
AR(4) process if they existed in the data. This process, which is a special case of an 
AR(4) process shall be specifically referred to as an AR(1,4) process. In addition, if 
the AR(1,4) process exists then there must be a way to transform the data and 
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calculate estimates for the parameters of the process. This proposed procedure for the 
AR(1,4) process would prevent oversights of existing processes and would save time 
expended in having to correct for the AR(1) process and then the AR(4) process. 
Most importantly though it would insure that if the AR(1,4) process was present in the 
data, it would be corrected. In this way any peculiarities in the data which might 
prevent detection and correction of either the AR(1) or AR(4) process would be 
avoided. 

As previously discussed, with an autoregressive disturbance present in the data, 
the OLS assumption of independent error terms is violated. This is related to problems 
in the error covariance matrix which is denoted as 

E(ee') = (eqn 1.14) 

where e is a vector of regression residuals. The problem in the error covariance matrix 
is that there exist off diagonal terms not equal to zero. This condition can indicate 
that some form of autocorrelation exists. An additional problem in the error 
covariance matrix is that the diagonal elements are not equal. This condition indicates 
that the error terms are not distributed with a constant variance. Errors distributed in 
this manner cause estimates of regression coefficients to be inefficient though unbiased. 
If the error covariance matrix equals 

(eqn 1.15) 



where I is an identity matrix, the regression coefficient P is correctly defined by 



The general procedure [Ref. 6: p. 440] used to attain an efficient estimate for P is: 

• Find a matrix P such that P'P = 

• Using the matrix P, transform the original data set where 



instead of 




(eqn 1.16) 



A 



(X'v^XrX'y^Y. 



(eqn 1.17) 



A 
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(eqn 1.18) 



Y* = PY 
and 

X* = PX. (eqn 1.19) 



• Perform GLS on the transformed model 

Y* = X*p + e* (eqn 1.20) 

where 

e* = Pe (eqn 1.21) 

to estimate the regression coefficients where 
/\ * * 1 * 

p = (X 'X )-‘X Y . (eqn 1.22) 



C. RESEARCH QUESTIONS 

Following the procedure as outlined above, this paper will address the AR(1,4) 
process. Specifically, this thesis will attempt to develop: 

• A statistic capable of detecting the AR(1,4) disturbance. 

• A P matrix for a transformation to account for the AR(1,4) disturbance. 

• A method to estimate parameters required for the transformation. 

D. RESEARCH METHODOLOGY 

In reaching the first objective, the statistic will be calculated by estimating the 
distribution of the ratio of two quadratic forms in normal variables. This method is 
based on ImhofF s [Ref. 7] technique and is outlined by Koerts and Abrahamse [Ref. 8]. 
These values shall be calculated for an extension of Schmidt's statistic [Ref. 9] to the 
AR(1,4) process. This statistic was developed in a manner similar to Durbin and 
Watson's and is essentially based on their work. The power of this statistic to detect 
the AR(1,4) process will be compared to the current sequential method of testing. 
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To account for the autoregressive disturbance, it is required that there be a P 
matrix capable of transforming the data with the property that 



The inverse of the y matrix for the fourth order process was derived by Siddiqui 
[Ref. 10]. With the form of the y inverse matrix known, a P matrix with the above 
stated property can be found using a method proposed by Beach and MacKinnon 
[Ref. 11]. Both the y inverse matrix and P matrix are expressed in terms of 0's or 
autoregressive parameters from the transfer function 



Utilizing the Yule-Walker equations [Ref. 12: p. 55], the autoregressive parameters 
from this transfer function were expressed in terms of the autocorrelation coefficients 
or p's. Once the 0's were determined in terms of the p's, substitutions were made into 
the derived P matrix. The resulting P matrix can be used to transform an AR(4) 
process as represented by the error term above. Since this error term is not 
representative of the AR(1,4) model, corrections were made to the P matrix by setting 
appropriate values of p equal to zero. For Schmidt's statistic, which will be introduced 
in Chapter II, p-> and p-j were set equal to zero. The resulting form of the P matrix is 
that required to transform data with an AR(1,4) disturbance. 

With the P matrix defined, the next task is to find estimates for pj and p^. Since 
the method outlined above is generalized least squares, one approach for estimating the 
values of pj and p^ is by the sample correlation coefficients 



where e t is a vector of least squares residuals. Since estimators derived in this manner 
are subject to considerable small-sample bias, the approach which will be developed 
will be to estimate pj and p^ by direct utilization of the test statistic. This method is 
analogous to Durbin and Watson's method where 



P'P = y' 1 . 



(eqn 1.23) 




+ n t 



(eqn 1.24) 



r s = I VtVl e t 2 s=1 > 4 



(eqn 1.25) 



Pj = 1 - .5dj 



(eqn 1.26) 
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and dj is the Durbin-VVatson statistic for an AR(1) process. 

After developing the above outlined procedure and generating the bounds for 
Schmidt's statistic a major question which must be answered is: How well do the 
procedures work in comparison to each other? In order to make the comparison, 
simulations will be used which utilize least squares residuals generated for various 
values of p j and p_| using a formula for the error of 

c t = Pl c t-1 + P4 e t-4 + n t ( ec l n *- 27 ) 

where T| t is distributed normal with zero mean and specified variance T . With 
residuals calculated as such, the Schmidt, Durbin-Watson and Wallis statistics will be 
calculated. In addition to the Durbin-Watson statistic, after appropriately 

transforming the residuals, the Wallis statistic will be used in the two-step method to 
test for the presence of an AR(4) disturbance. For various sets of residuals generated 
with different variances for the normal error term, the power of the statistics to detect 
the AR(1,4) disturbance will be determined. Additionally, the number of observations, 
T, will be varied by varying the number of residual terms. Thus the power of the 
various statistics to detect the AR(1,4) process will be determined for different values 
of T, pj and p_j and x . Residuals used in the above calculations will also be used to 
calculate estimates of pj and p^. Utilizing mean square error, calculated from 
simulation data, the efficiency of derived formulas to estimate Pj and p^ will be 
determined. Specifically, the ability of both procedures to detect the AR(1,4) 
disturbance will be determined for a best case with T= 100 and x — .01 and a worst 
case where T = 20 and T 2 = 10. 

The ability of the entire procedure derived for Schmidt's statistic in estimating the 
regression coefficients will be determined for a best and worst case scenario. The best 
case will be for t 2 = .01 and T = 100, while the worst case will be for T 2 = 10 and T 
= 20. Using the regression model 

Y t = px t + £ t (eqn 1.28) 

where 

c t = Pl £ t-1 + Pt-4 + V ( e fi n 1 - 29 > 
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for X t = .1 and P = 1 the Y t s are generated for various values of pj and p_j. Once 
the Y t 's are determined, an estimate of P will be calculated. Mean square error will be 
used as a measure to compare estimates of P with their known value of one. This 
measure will indicate the efficiency of regression parameters estimated by the procedure 
as based on Schmidt's statistic. 

E. THESIS ORGANIZATION 

The thesis is organized into four chapters. Chapter II will introduce background 
theory' utilized to derive the AR(1,4) procedure. In particular, the Durbin- Watson test 
for the AR(1) process will be addressed along with the Wallis and Schmidt tests for 
special cases of the AR(4) process. In addition to test procedures, the transformation 
for autocorrelated data will be discussed. 

Chapter III will initially develop the test for the AR(1,4) process. Specifically, 
upper and lower bounds will be developed for the Boger and Schmidt statistics. After 
deriving the test, the data transformation for the AR(1,4) process will be addressed. 
Next, results will be presented concerning the effectiveness of the derived results for the 
AR(1,4) process. Finally, conclusions and areas of future study will be discussed. 
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II. BACKGROUND THEORY 



A. GENERAL 

In this section certain concepts underlying all time series processes will be 
defined. These definitions will lay the groundwork which will lead into discussion of 
specific test procedures. In particular, the Durbin-Watson, Wallis and Schmidt tests 
shall be discussed. Discussion of these procedures will act as a prelude to the next 
chapter's development of a test and correction procedure utilizing an extension of 
Schmidt's statistic. 

A sample of error residuals ej,....ey, where the index denotes a point in time, is 
called a time series. Each observation in this sample is a realization of a random 
variable and as a result this sequence is considered to be a discrete stochastic process. 
The sequence is considered discrete since observations are made at a fixed interval of 
time. Additionally, ej,....ey is considered as a finite subsequence of an infinite series 
whose joint distribution is defined by these finite subsequences. 

Only stationary stochastic processes will be considered in this paper. A process is 
considered stationary if the joint probability distribution of the residuals ej,....ey is 
unaffected by shifting the time origin forward or backward by k time units. That is, 
ej,....ey and e^ + i>— <e k + T’ are identically distributed [Ref. 12). 

The general linear model with autoregressive disturbance of order p is defined as 

Y t = X t 'p + c t (eqn 2.1) 

with 

c t = 6 l c t-l + - + 6 p c t-p + ^t ( ec l n 2 - 2 ) 

where Y t is the value of the dependent variable at time t, X t ' is a (1 x K) vector 
representing an observation on K non-stochastic variables at time t and P is a (K x 1) 
vector of coefficients to be estimated. Equation 2.2 implies that error disturbances in 
the current time period depend on those in previous time periods and another error 
term which is independent with zero mean and constant variance. 
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Autocovariance defines the linear relationship between members of a stochastic 
sequence. In particular, the autocovariance, y^, between c t and + ^ is given by: 

Yk = E (( £ t' E ^ c t)X c t+k’ E ( £ t + k))) k = 0,1,2,... (eqn 2.3) 

where E(£ t ) is the mean of the stochastic process, which is zero in our application. For 
a stationary 7 process the covariance between £ t and £ t +k depends on k and not on the 
particular point in time. The above formula completely describes the autocovariance 
structure of the stochastic sequence as defined above [Ref. 13: p.226]. 

Autocovariance is dependent on the unit of measurement for the underlying 
variable. To account for this, the y^'s are normalized by dividing by y*, which is the 
variance of e v Dividing y^ by */q results in the autocorrelation function of £ t : 

Pk = Yk/Yo> k= 0,1,2,.. (eqn 2.4) 

From the above definition of the autocorrelation function it is obvious that 

p 0 = 1. (eqn 2.5) 

In order to gain estimates for the 0's or autoregressive parameters in terms of the 
autocorrelation coefficients or p's, the Yule- Walker equation is used [Ref. 12: p. 55). 
This equation for the above general linear model is derived by left multiplying the 
equation 

£t = e i £t.i + ... + e p £t.p+nt (eqn 2.6) 

by £,_k> an d then taking expectations and dividing by the variance Jq of £ r This 
results in the Yule- Walker equation 

Pk = e iPk-l + - + 0 pPk-p fork* 1,2 (eqn 2.7) 

By substituting k= 1,2, ...p into the above equation, a set of linear equations for 
0j,02- -0^ in terms of Pj^-.-Pp are obtained. By substituting estimates for P|.p->...Pp 
into the above equations, the autoregressive parameters are estimated. Various other 
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methods exist to estimate the p's, and some will be discussed later. This result is 
important in the specification of the autoregressive process and will help in the 
development of a transformation matrix for the autocorrelated data [Ref. 12]. 

B. DURBIN- WATSON TEST 

This section summarizes the most important features of the theory for the 
Durbin-Watson test (Ref. 2,14,15]. This theory is the basis for later generalizations of 
tests for autoregressive disturbances. 

The Durbin-Watson test is based on the statistic 



d = I( e r e t-i) 2 /I c t 2= e ' Ae / e ' e 

<*y 

where e = y - XP is a vector of least squares residuals and A equals: 



(eqn 2.8) 



1 -1 0 
-1 2 -1 
0-12 



0 0 0 

0 0 0 

-10 0 



0 0 0 
0 0 0 
0 0 0 



0* 

0 

0 



0 0 0 -1 2 -1 
0 0 0 0 -1 1 

Ori 



For certain cases where regression vectors are eigenvectors of the matrix A 
occurring in the residuals distribution, the statistic 



e'Ae / e’e (eqn 2.9) 

provides a test which is uniformly most powerful against certain alternative hypotheses 
[Ref. 16: p. 88]. 

The general linear model was previously defined as 
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Y = XP + £, 



(eqn 2.10) 



where Y is a (T x 1) matrix, X is a (T x K) matrix of T observations on K variables 
and c is a (T x 1) matrix of errors. The least squares estimate of p is P which is given 
by: 

(X'X^X'y. (eqn 2.11) 

Also, e, the vector of residuals from the regression, is defined as 

e = y - X$ = My (eqn 2.12) 

where I is an identity matrix of order T, and 

M = I - X(X'X)‘ 1 X'. (eqn 2.13) 

It can be verified that the matrix M is idempotent. Substituting XP + e for y into the 

formula for e leads to the result 

e = Me. (eqn 2.14) 

This result implies that 

d = e'MAMe / e'Me. (eqn 2.15) 

The distribution of the above statistic is dependent on the (T - K) non-zero eigenvalues 
of AM which are denoted by S-. Values for Sj are dependent on M through the matrix 
X. 

The matrix A is a (T x T) symmetric matrix with (T - 1) positive eigenvalues 
which for T greater than or equal to three lie in the closed interval from zero to four. 
The eigenvalues of A will be denoted by Xj and are defined as 

X- = 2(1 - cos(rcj/T)) j = 1 T-l. (eqn 2.16) 
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Durbin and Watson showed that if the e t are assumed to be normally distributed, 
there exists an orthogonal linear transformation of e to v such that 

d = Z 6 i v i 2 / E v i 2 < (eqn 2.17) 

tz/ i»i 

where Vj are independent and identically distributed normal variables with zero mean 
and variance <x 2 . The distribution of d thus depends on the which are the 
eigenvalues of M [Ref. 2: p. 412]. 

Using the above result to obtain the distribution of d is tedious, since the 8j's 
depend on the X matrix and the statistic would have to be calculated for each new 
matrix. Consequently, if it is desired to have a test which is not restricted to a 
particular X matrix, the distribution would have to be determined for all X matrices. 
This task is impossible since there are an infinite number of matrices. 

Durbin and Watson avoided this problem by constructing a bounds test where 
the upper and lower bounds, d U and d^, respectively, are independent of the particular 
X matrix. They were able to accomplish this by determining inequalities on the in 
terms of the eigenvalues Xj of the matrix A, where the 8-'s and Xj's have been arranged 
in increasing order. This result implies that 

d L <d<d y (eqn 2.18) 

where 

d U = X X i+ K v i 2 / I* v i 2 > ( ec l n 2 - 19 > 

(•/ Ut 

and 

T-K . r--K 

d L = IVi / X V i • ^ ec l n 2 - 20 > 

/X / * 

Durbin and Watson were unable to find the exact distribution of these bounds but 
approximated the distribution with the first four terms of a series expansion in terms of 
Jacobi polynomials using a Beta distribution as a weight function. 

Using the statistic 

t - r 0 

X( e t - e t -i > / S e t ’ ( ec i n 2 - 21 ) 

-t-/ 
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and asymptotic results it can be shown that 



d = 2( 1 - p ), 



(eqn 2.22) 



where p is an estimate of the autocorrelation coefficient. This result shows that when 
p equals zero, indicating no autocorrelation, d equals two. In addition, when p equals 
one indicating positive autocorrelation, d equals zero and when p equals negative one, 
d equals four. In testing the null hypothesis of no autocorrelation against the 
alternative hypothesis of positive autocorrelation, the null hypothesis would be rejected 
if the calculated value for d is less than the lower bound calculated and accepted if d is 
greater than the upper bound. If the calculated value of d falls between the upper and 
lower bounds the test is considered inconclusive and there is insufficient evidence to 
accept or reject the null hypothesis. 

Since Durbin and Watson first calculated upper and lower bounds for their 
statistic using an approximate distribution, a method has been developed by Imhoff 
which allows the calculation of an exact distribution for a quadratic form in normal 
variables [Ref. 7: p. 419|. To briefly indicate how the Imhoff distribution is utilized in 
the calculation of the statistic's bounds, only the lower bound will be considered. The 
calculation of the upper bound can be done in an analogous manner. In order to 
calculate the lower bound, the eigenvalues of the A matrix, X-, are calculated. For the 
A matrix, there will be (T - 1) non-zero eigenvalues where T is the order of the A 
matrix and the one zero term is the eigenvalue which corresponds to the constant term 
of the regression. From these eigenvalues, the (T - K) smallest are selected to compute 
the statistic's lower bound [Ref. 8: p.70[. The value for K is equal to the number of 
regressors including a constant term. The basis for the above procedure is the result 
proved by Durbin and Watson which determined inequalities of the eigenvalues Sj of 
the matrix AM in terms of the eigenvalues X- of the A matrix or 



i = 1,2,...T-K , 



(eqn 2.23) 



where 



K' = K - 1 



(eqn 2.24) 
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and K' is the number of regressors not including a constant term. Analogously, to 
calculate the upper bound the (T - K) largest eigenvalues are selected. Adapting the 



Imhoff distribution to this problem leads to 

F(d L ) = 1/2 - 1/rcj (( sin <p(u)) / (uto(u))) du 

0 

where 


(eqn 2.25) 


7-K . 

<p(u) = .5 J tan'Tuju) 
/«/ 


(eqn 2.26) 


w(u) = ri(i + u:^u^)^ 
/!/ 

and 


(eqn 2.27) 


uj = \ - d L (i = 1 T-K) 


(eqn 2.28) 



Since the Xj's are known, a value for d^ is selected and the uj's are calculated. 
Using the Uj's, the integral is then evaluated numerically. Numerical integration of this 
infinite integral is possible since Imhoff has proved that the integral's limit as u 
approaches zero is 



T-K 

l/2£ up 
/*/ 


(eqn 2.29) 



Since the integral's range is infinite, there will be a truncation error associated with this 
integration. Imhoff has shown that the truncation error caused by integrating over the 
finite range from 0 to u inclusive, can be held to any arbitrary level p by taking 



u = (((T - K)/2) rcpnVjl 1 /^ )g 
/-/ 


(eqn 2.30) 



where g = 2/(T-K). Additionally, the numerical integration will be subject to error 
related to the method of integration. This error is controlled by setting the error 
tolerance within the program used to integrate the function. 
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C. EXTENSIONS TO HIGHER ORDER PROCESSES 

Durbin and Watson's original work has been extended to higher order 
autoregressive processes. Specifically, Wallis has extended it to a special case of an 
AR(4) process and Schmidt has extended it to a second order process. For the 
regression model 

Y = XP + £, (eqn 2.31) 

with an error specification 

c t = P £ t-4 + H t » (eqn 2 - 32 ) 

Wallis generated bounds for the statistic 

r r 

d 4 = X( e t * e t-4^ /X e t = e A 4 e / e ' e - ( ec l n 2 - 33 ) 

Wallis' derivation of this result is analogous to Durbin and Watson's. Using 
ImhofTs distribution, Wallis was able to calculate bounds for this special case of the 
AR(4) process [Ref. 5: p. 617], Vinod, in a manner similar to Wallis, further extended 
these results to include bounds for the same special cases of the AR(2) and AR(3) 
processes [Ref. 17]. 

In order to account for an AR(1) or AR(4) disturbance, the data must be 
transformed. In this regard, for the AR(1) process the dependent variables are 



transformed as 

Y t * = YjO-pj 2 ) 1 / 2 , t=l (eqn 2.34) 

and 

Y* = Y t - Pl Y M , t= 2,3 T. (eqn 2.35) 

Similarly, the independent variables are transformed as: 

X t * = X t (l - Pj 2 ) 1 / 2 , t=l (eqn 2.36) 
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and 



X t * = X t - pjX t _!, t = 2,3, T. (eqn 2.37) 

In order to utilize this transformation, an estimate for pj is required. Several 
methods are available to estimate this parameter. The method of choice due to its 
favorable properties and ease of computation [Ref. 3] is 

Pj = 1 - .5dj. (eqn 2.38) 

For the special case of the AR(4) process the transformation is similar. 
Specifically, the dependent variables are transformed as: 



Y t * = Y t (l - p/) 1 / 2 t= 1,2,3, 4 (eqn 2.39) 

and 

Y t * = Y t - p 4 Y t _ 4 t = 5 T. (eqn 2.40) 

The independent variables, similarly are transformed as: 

X t * = X t (l - p/) 1 / 2 t= 1,2, 3, 4 (eqn 2.41) 

and 

X t * = x t - p 4 X t . 4 t = 5, T. (eqn 2.42) 



As in the AR(1) case, an estimate of p 4 is required. The formula for this estimate is 
the same as before except that the Wallis statistic is used in lieu of the Durbin-Watson. 
Once either of these transformations is utilized after the detection of its respective 
disturbance, estimates of regression parameters can be expected to be efficient and 
unbiased. 

Schmidt [Ref. 9] generalized Durbin and Watson's preliminary work to the case 
of a second order autoregressive process whose disturbance is of the form: 
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£ t = Pl £ t-1 + f>2 £ t-2 + H t - 



(eqn 2.43) 



The statistic which Schmidt derived to detect this disturbance was: 



d 2 = d i + 6 2 

where 



= 2(« t 

i*L 



- e 



t- 




2 

t 



and 



s 2 = 2( e t - e t-2> 2 / S e t 2 ' 

In matrix notation this statistic is represented as 



(eqn 2.44) 



(eqn 2.45) 



(eqn 2.46) 



d 2 = e'(Aj + A 2 )e / e'e (eqn 2.47) 

where Aj is the matrix as previously defined for the Durbin-Watson statistic and A 2 is 
defined as 



1 0 -1 
0 1 0 
-10 2 
0-10 



0 0 0 

-10 0 
0-10 
2 0-1 



0 0 0 
0 0 0 
0 0 0 
0 0 0 



0 

0 

0 

0 



0 -1 0 2 0 -1 

0 0 10 10 
0 0 0 1 0 1 
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Following results proved by Durbin and Watson, Schmidt generated the statistic, d 2 , 
utilizing ImhofFs distribution and the eigenvalues from the matrix (Aj + A 2 ). In a 
manner similar to previous developments, he further showed that 



<^2 “ 2(2 - Pj - P2)- 



(eqn 2.48) 



D. SEQUENTIAL TESTING 

In development of a generalization of the Durbin-Watson statistic, Vinod, 
derived a scheme of sequential tests for the presence of autocorrelation [Ref. 17]. In 
this procedure a sequence of tests is used to determine if AR processes of sequentially 
increasing order exist. As an example, the null hypotheses for the four tests required 
to detect a fourth order process would be: 

• H 0 : pj = 0, 

• Hq : p 2 = 0 given p } = 0, 

• Hq : p 3 = 0 given Pj = P 2 = 0, 

• H 0 : p 4 = 0 given Pj = p 2 = P3 = 0. 

Thus in order to test for an autoregressive disturbance, the null hypotheses for 
sequentially greater orders are tested until one of the null hypotheses is rejected. 
Specifically, at the step, j > 1, the hypothesis, Hq : Pj = 0 given that Pj = p 2 = 
... = Pj_j = 0 is tested against the two-sided alternative H a : Pj > 0 and Pj < 0. 
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III. THE AR(1,4) PROCESS 



A. INTRODUCTORY ANALYSIS 

In this section stationary conditions and the development of the error term for 
the AR(1,4) process will be presented. If stationary conditions do not hold the 
AR(1,4) model is invalid and a model which considers a dynamic specification must be 
pursued. In particular, stationary conditions must be determined for the regression 
model 



Y t =X t 'P + c t , (eqn 3.1) 

whose error term is given by: 

® l^t- 1 ^4^t-4 ^t’ (eqn 3.2) 

As discussed in previous chapters, the Tj t 's are assumed to be distributed normally with 
zero mean and constant variance. 

Wise [Ref. 18] derived a method to determine the stationary conditions for a 
stochastic process of the autoregressive and moving average types. Following his 
procedure, it is possible to show that the stationarity conditions for a fourth order 
autoregressive process are 

• 2 + pj - pj + 2p 4 > 0, 

• 3 + p 2 - 3p 4 > 0, 

• 2 - Oj + pj + 2p 4 > 0, 

• 1 - Pi - P2 - P3 - P4 > °- 

• 1 + Pi - P 2 + P3 - P4 > °- 

Since P 2 and p^ are assumed to be zero in the AR(1,4) process, these stationarity 
conditions reduce to 

• 2 + Pj + 2p 4 > 0, 

• 3(1 - p 4 ) > 0, 

• 2 - pj + 2p 4 > 0, 

• 1 - Pi - P 4 > 0, 

• 1 + Pi - Pq > 0. 
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The error process for a general fourth order process is defined by the relation: 



E t = 0 l E t-l + 0 2 E t-2 + 0 3 E t-3 + 0 4 E t-4- ( e< I n 3l3 ) 

If the relationship is multiplied through by and expectations are taken, the 
difference equation 

Yk = e iYk-l + 0 2Yk-2 + - + 0 4Yk-4 k>0 (eqn 3.4) 

results, where y was previously defined as the autocovariance. If this equation is 
divided through by Jq, the result is 



Pk -0 lPk-l +0 2Pk-2 + " +0 4Pk-4’ ( ec l n 3 -^) 

which is the general form of the Yule-Walker equation and defines the autoregressive 
parameters, 0's, in terms of the autocorrelations, p's. If k = 1 , 2, 3, 4, is substituted 
into this equation a set of linear equations is obtained for 0 j, © 2 , ©3 and 0 ^ in terms of 
Pj, p->, pj and p^. Specifically, these equations are 



Pj = 0j +02Pj +02P2 + 0 4 P3, (eqn 3.6) 

p 2 = e iPi+02 + e 3Pi+ e 4 P 2 . (eqn 3.7) 

P3 = 0 lP2 + 0 2 p l +0 3 + 0 4Pl’ (eqn 3.8) 

and 

P4 = 0 lP3 + 0 2P2 + 0 3Pl + 0 4 - (eqn 3.9) 



Since the error process for the AR(1,4) process is specified as 



E t = e i E M +0 4 E t-4 + n t > 



(eqn 3.10) 
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where ©2 and 63 are assumed equal to zero, the above equations can be reduced to 



Pj = 0j + 04 P3, (eqn3.11) 

P2 = 6 1 P 1 + 6qP2’ (eqn 3.12) 

P3 = © 1P2 + ®qPl> (eqn 3.13) 

and 

P4 = 0 1P3 + 64- (eqn 3.14) 

Further since it is assumed that P[ and P3 are not present in the disturbance, the above 
equations can be finally reduced to 

Pj = 0j (eqn 3.15) 

and 

P4 = 64 - (eqn 3.16) 

These final equations allow the AR(1,4) error term to be represented as: 

"b Pq£j _4 "l" (eqn 3,17) 

B. BOGER'S STATISTIC 

Boger proposed a statistic which was a variation of the Durbin-Watson statistic 
to test for the AR(1,4) disturbance. This statistic is of the form 

d 1 4 = E (e t -e t . j -e t . 4 ) 2 / £e t 2 . (eqn 3. 1 8) 

*•5 t 1 */ 
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In matrix notation this statistic is denoted as 



d 1<4 = e'A 1)4 e / e'e 
where A j 4 equals 



(eqn 3.19) 



1 0 0 

0 1 0 

0 0 1 

1 0 0 

-1 1 0 

0 -1 1 



1 -1 0 

0 1 -1 

0 0 1 

1 1 0 

-1 2 -1 

0-12 



0 0 0 

0 0 0 

1 0 0 

0-10 
0 0 0 

-10 0 



0 



- 110 - 12-1 
0-1 0 0-1 1 



Bounds for Boger's statistic could not be determined because its A matrix as 
defined above did not satisfy properties as established by Durbin and Watson. In 
particular, Durbin and Watson [Ref. 2] showed that the matrix product MA, where A 
is any real symmetric matrix has (T-K) real positive eigenvalues of which K are zero 
eigenvalues. T and K are the dimensions of the matrix X. Specifically, Durbin and 
Watson showed for their statistic whose A matrix equals 



1 


-1 


0 


0 


0 






0 


-1 


2 


-1 


0 


0 








0 


-1 


2 


-1 


0 








0 
















0 








. 








0 


. 






. 








0 










2 


-1 


0 


0 










-1 


2 


-1 


0 


• 


• 






0 


-1 


1 
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that there were T-l positive eigenvalues. As previously discussed, they also showed 
that if the positive eigenvalues from both MA and A are placed in increasing order, 
bounds for statistics of the form 



d = e'Ae / e'e, 



(eqn 3.20) 



can be calculated. 

Where Boger's statistic violated these properties was that for all values of K or 
the number of independent variables, there were not sufficient eigenvalues in the A 
matrix to bound those eigenvalues in the matrix MA. Specifically, the A matrix for 
Boger's statistic had one zero and three negative eigenvalues where Durbin and 
Watson's had one zero eigenvalue. Since the matrix MA had T-K positive eigenvalues 
and Boger's A matrix only had T-4 positive eigenvalues instead of T-l, when the 
number of independent variables was less than or equal to three, there were not 
sufficient eigenvalues from Boger's A matrix to bound those from the matrix MA. 
Since this violated Durbin and Watson's theory, upper and lower bounds could not be 
calculated for Boger's statistic, when the number of independent variables was less than 
three. 

C. SCHMIDT'S STATISTIC EXTENDED TO THE AR(1,4) PROCESS 

Since it w r as desirable to have a statistic which was not limited by the number of 
independent variables, a variation of Schmidt's statistic to fit the AR(1,4) process w'as 
investigated. Schmidt was able to show that the statistic d^ where 

d 2 = dj + <>2 (eqn 3.21) 



provides a test to detect a second order autoregressive disturbance. 

Using Schmidt's research this paper proposes that it is possible to detect the 
AR(1,4) process by way of the statistic dj 4 where: 



and 




(eqn 3.22) 



d 1 ,4 = d ! + d 4 



(eqn 3.23) 



31 



or 



In matrix notation this statistic is given as 




(eqn 3.24) 



e'Aj 4 e / e'e 



(eqn 3.25) 



where 



A l,4 - Aj + A 4 



(eqn 3.26) 



or A j 4 equals: 

"2 -1 0 0 -1 0 0 0 0 0 

-1 3 -1 0 0 -1 0 0 0 0 

0 -1 3 -1 0 0 -1 0 0. 

0 0 -1 3 -1 0 0 -1 0 0 

-1 0 0 -1 4 -1 0 0 -1 0 



Results established by ImhofF and Durbin and Watson were used to determine 
upper and lower bounds for the statistic dj 4 . Koerts and Abrahamse [Ref. 8] have 
provided Fortran 66 code which allows the determination of the distribution of a 
quadratic form in normal variables by way of ImhofFs distribution. This code was 
updated to Fortran 77 standards and modified to allow the calculation of the dj 4 
statistic. The algorithm utilizes ImhofFs distribution to determine five percent 
significance points of the statistic by successively halving the range of the function 
until 



-1 0 0-1 3-1 



-10 0-12 



F(d L ) = .05, 



(eqn 3.27) 



32 



where F(dj^) is the Imhoff distribution function as defined in the previous chapter. 
Truncation and integration errors are controlled by input parameters. It was found 
that if these errors were set to less than .0001, the program required inordinate 
amounts of CPU time. At a .0001 level of accuracy, the five percent significance points 
utilized S2500 of computer resources. At a .00001 level of accuracy, approximately a 
six fold increase in resources can be expected. Five percent significance points for the 
statistic dj ^ are listed in Table I. In Table I, K is the number of regressors including 
the intercept term and T is the number of observations. 

Currently, if an AR(1,4) process is suspected to exist in the data a two-step 
procedure is generally used. Specifically, this involves first testing and correcting for 
the AR(1) component of the disturbance and then doing the same for the AR(4) 
component. The AR(1,4) procedure proposes to test and correct the AR(1,4) 
disturbance in one step where both components are simultaneously corrected. Once 
autocorrelation has been determined to exist within the data for both procedures, 
values for the autocorrelation coefficients Pj and p 4 must be determined. In the 
two-step procedure, Pj is calculated by way of the formula 



where dj is the Durbin- Watson statistic. After an appropriate transformation, the data 
are tested in a sequential manner as discussed in the previous chapter until a higher 
order disturbance is detected. If a fourth order disturbance is detected, p 4 is estimated 
by the formula 



where d 4 is the statistic tabulated by Wallis. As can be seen by this procedure, the 
calculation of values for p j and p 4 are in a sense conditional upon each other. That is 
since pj's value is used as part of the original data transformation and p 4 is determined 
by way of the Wallis statistic which is calculated from the transformed data, the value 
of p 4 depends on the value of pj. 

There are two important assumptions concerning the two-step procedure as it 
relates to the AR(l,4) disturbance. The first of these is that the Durbin-Watson 
statistic will be used to detect the AR(1,4) process. That is, that an individual who is 



Pj — 1 - .5dj, 



(eqn 3.28) 



p 4 = 1 - .5d 4 



(eqn 3.29) 
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TABLE I 

TABLE OF BOUNDS FOR THE SCHMIDT STATISTIC 



K -- 



4 



T 


L 


U 


L 


U 


L 


U 


15 


2.38 


2.84 


2.19 


3.05 


2.03 


3.29 


16 


2.43 


2.87 


2.26 


3.06 


2.10 


3.30 


17 


2.48 


2.91 


2.31 


3.10 


2.17 


3.30 


18 


2.53 


2.94 


2.36 


3.13 


2.23 


3.30 


19 


2.57 


2.97 


2.42 


3.15 


2.27 


3.32 


20 


2.61 


3.00 


2.47 


3.16 


2.33 


3.34 


21 


2.65 


3.02 


2.51 


3.18 


2.37 


3.35 


22 


2.68 


3.03 


2.55 


3.20 


2.42 


3.35 


23 


2.72 


3.05 


2.58 


3.22 


2.46 


3.36 


24 


2.75 


3.08 


2.62 


3.24 


2.50 


3.37 


25 


2.78 


3.09 


2.65 


3.26 


2.54 


3.38 


26 


2.80 


3.11 


2.69 


3.26 


2.57 


3.39 


27 


2.82 


3.12 


2.71 


3.28 


2.60 


3.40 


28 


2.85 


3.14 


2.74 


3.29 


2.63 


3.41 


29 


2.87 


3.15 


2.77 


3.30 


2.66 


3.42 


30 


2.90 


3.17 


2.80 


3.31 


2.69 


3.43 


31 


2.92 


3.17 


2 . 82 ' 


3.32 


2.72 


3.44 


32 


2.94 


3.19 


2.84 


3.33 


2.74 


3.44 


33 


2.96 


3.20 


2.86 


3.33 


2.77 


3.45 


34 


2.98 


3.21 


2.88 


3.34 


2.79 


3.46 


35 


2.99 


3.22 


2.90 


3.35 


2.81 


3.46 


36 


3.01 


3.23 


2.92 


3.36 


2.83 


3.47 


37 


3.03 


3.24 


2.94 


3.36 


2.85 


3.48 


38 


3.04 


3.25 


2.96 


3.37 


2.87 


3.48 


39 


3.05 


3.26 


2.97 


3.38 


2.89 


3.49 


40 


3.07 


3.27 


2.99 


3.38 


2.91 


3.49 


45 


3.13 


3.31 


3.05 


3.41 


2.98 


3.50 


50 


3.18 


3.34 


3.12 


3.43 


3.05 


3.52 


55 


3.23 


3.37 


3.16 


3.46 


3.10 


3.54 


60 


3.26 


3.40 


3.21 


3.47 


3.15 


3.55 


65 


3.30 


3.42 


3.24 


3.49 


3.19 


3.56 


70 


3.33 


3.44 


3.28 


3.51 


3.23 


3.57 


75 


3.35 


3.46 


3.31 


3.52 


3.26 


3.58 


80 


3.38 


3.47 


3.34 


3.53 


3.30 


3.59 


85 


3.40 


3.49 


3.36 


3.54 


3.32 


3.59 


90 


3.42 


3.50 


3.38 


3.55 


3.34 


3.60 


95 


3.43 


3.51 


3.40 


3.56 


3.36 


3.61 


100 


3.45 


3.53 


3.42 


3.57 


3.38 


3.62 
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TABLE I 

TABLE OF BOUNDS FOR THE SCHMIDT STATISTIC (CONT'D.) 



K-- 5 6 



r 


L 


U 


L 


U 


15 


1.87 


3.47 


1.73 


3.74 


16 


1.95 


3.47 


1.82 


3.73 


17 


2.01 


3.47 


1.88 


3.72 


18 


2.07 


3.48 


1.95 


3.71 


19 


2.13 


3.49 


2.01 


3.70 


20 


2.19 


3.49 


2.07 


3.68 


21 


2.24 


3.50 


2.12 


3.67 


22 


2.29 


3.51 


2.17 


3.66 


23 


2.33 


3.52 


2.22 


3.66 


24 


2.38 


3.52 


2.27 


3.65 


25 


2.42 


3.52 


2.31 


3.64 


26 


2.45 


3.53 


2.35 


3.64 


27 


2.49 


3.53 


2.39 


3.65 


28 


2.52 


3.54 


2.42 


3.65 


29 


2.56 


3.54 


2.45 


3.65 


30 


2.59 


3.54 


2.49 


3.65 


31 


2.61 


3.55 . 


2.52 


3.65 


32 


2.65 


3.55 


2.55 


3.65 


33 


2.67 


3.55 


2.58 


3.65 


34 


2.70 


3.56 


2.61 


3.65 


35 


2.72 


3.56 


2.63 


3.66 


36 


2.75 


3.57 


2.66 


3.66 


37 


2.77 


3.57 


2.68 


3.66 


38 


2.79 


3.57 


2.70 


3.66 


39 


2.80 


3.57 


2.72 


3.66 


40 


2.83 


3.58 


2.74 


3.66 


45 


2.91 


3.59 


2.84 


3.66 


50 


2.98 


3.60 


2.92 


3.67 


55 


3.05 


3.62 


2.99 


3.68 


60 


3.10 


3.62 


3.04 


3.69 


65 


3.15 


3.63 


3.09 


3.69 


70 


3.18 


3.63 


3.13 


3.69 


75 


3.22 


3.64 


3.18 


3.69 


80 


3.25 


3.64 


3.21 


3.70 


85 


3.28 


3.65 


3.24 


3.70 


90 


3.30 


3.66 


3.27 


3.71 


95 


3.33 


3.66 


3.29 


3.71 


100 


3.35 


3.66 


3.32 


3.71 
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testing for the AR(1,4) process will not calculate the Wallis statistic if the 
Durbin-Watson statistic fails to detect. This assumption should be fairly accurate in 
practice as an individual would most likely decide that autocorrelation is not present if 
the Durbin-Watson statistic indicated such. The second assumption which is closely 
related to the first, concerns the order of calculation for pj and p 4 . Specifically, as 
previously shown, the formula derived from the Schmidt statistic contains a Pj and p 4 
term. Since this is the case, different results might have been obtained if p 4 had been 
calculated by way of the Wallis statistic formula first, after which Pj would have been 
calculated by the Schmidt statistic formula. This possibility w'as not investigated in 
this paper, but Pj was calculated by way of Durbin and Watson's statistic after which 
p 4 was determined from the Schmidt formula. 

For the AR(1,4) procedure values for Pj and p 4 wall be calculated in a similar 
manner. Specifically, it is proposed that Pj will be calculated by the formula 

Pj = 1 - .5dj (eqn 3.30) 

and p 4 by the formula 

p 4 = 2 - pj - .5dj 4 . (eqn 3.31) 

The formula for p 4 is derived from the original definition of dj 4 or 



If the formula for dj 4 is factored and asymptotic arguments applied, after appropriate 
rearrangement of terms, the formula for p 4 is achieved. The validity of these formulas 
is argued in a manner similar to that for the two-step procedure. The difference is that 
values for pj and p 4 are determined in a single step. An advantage to this method is 
that a value for p 4 is calculated prior to the data transformation, thus its value will not 
be influenced by any peculiarities which may have occurred during the data 
transformation. 

Once estimates for pj and p 4 have been determined, a method to utilize them in 
correcting the AR(1,4) disturbance must now r be determined. As previously discussed, 
this requires that a matrix P be found such that P'P = ly'* where is the inverse of 




(eqn 3.32) 
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the sample covariance matrix. Once the P matrix is determined, the data can be 
transformed such that 



Y* = PY, (eqn 3.33) 

X* = PX, (eqn 3.34) 

e = Pe, (eqn 3.35) 

and least squares applied to the model 

Y = X P + e , (eqn 3.36) 

where the estimator for p is 

£ = (X'P'PX)' I X'P'PY. (eqn 3.37) 



If P is determined from the transformed data, it should be more efficient than estimates 
from ordinary least squares. 

Beach and MacKinnon [Ref. 11], in their paper on maximum likelihood 
estimation of a second order process, defined the P matrix as 



a 



b c 

-p 2 -pj 1 



O' 

0 

0 



Lo 



-p 2 -p l i 
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They proposed a secondary method to determine values for a, b and c where P'P is 
constrained to be proportional to the inverse of the variance-covariance matrix as 
given be Siddiqui [Ref. 10] and a, b and c are solved for in the implied restrictions. In 
a manner similar to that proposed by Beach and MacKinnon, this paper will solve for 
the P matrix of the fourth order process. 

From Siddiqui's paper on the inversion of the sample covariance matrix it was 
possible to determine that the inverse of the sample covariance matrix, V|/'*, equals 



i 



o 

0 






-*1 




W 



0 




0 



0 

0 




-% 



V* 

0 



2 2 



o 



0 



0 

0 

0 

0 

0 



0 



0 

0 

0 

0 

0 



0 

0 
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where a j = p j and a^ = p^. This matrix was obtained from the vj/“ ^ matrix for a fourth 
order process after setting the terms for a 2 and a^ equal to zero, where a^p-, and 

a 3 = f*3* 

The P matrix for the AR(1,4) process is defined such that P equals: 
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a 



. 0 

be 0 

def 0 

ghij 0 

-P 4 0 0 -Pi 1 0 

iP *P4 0 0 *Pj 1 

By following Siddiqui's method in the context as proposed by Beach and MacKinnon, 
10 equations in 10 unknown are defined 

a 2 + b 2 + d 2 -l- g 2 = 1 - p 4 2 , (eqn 3.38) 

be + de + gh = -p i , (eqn 3.39) 

fd + ig = 0, (eqn 3.40) 

gj = -Pjpq, (eqn 3.41) 

c 2 + e 2 + h 2 = 1 + pj 2 - p 4 2 , (eqn 3.42) 

fe + hi = -p j, (eqn 3.43) 
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jh = 0 , 



(eqn 3.44) 



f 2 + i 2 = 1 + p , 2 - p 4 2 , 


(eqn 3.45) 


ij = -Pi, 


(eqn 3.46) 


i 2 = 1 * P 4 2 


(eqn 3.47) 



Solving these equations for the 10 unknowns, the elements of the previously defined P 
matrix are determined where 



a = V (1 - P 4 2 - b 2 - d 2 - g 2 ), 


(eqn 3.48) 


b = (-pj - gh -de)/ c, 


(eqn 3.49) 


c = V (1 + pj 2 - p 4 2 - h 2 - e 2 ), 


(eqn 3.50) 


d = (-Pi 2 P 4 V(1 - P 4 2 ))/((l - P 4 2 ) V ((1 - P 4 2 ) 2 - Pj 2 P 4 2 )), 


(eqn 3.51) 


e = (-PjV (i - P 4 2 )) / (V ((1 - P 4 2 ) 2 - Pi 2 P 4 2 )), 


(eqn 3.52) 


f = V (((1 - P 4 2 ) 2 - Pj 2 P 4 2 ) / (1 - P 4 2 )), 


(eqn 3.53) 
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g = -P1P4 / J (1 - P4 2 ), 



(eqn 3.54) 



h = 0, 



(eqn 3.55) 



i = -Pi / V (1 - Pq 2 ), 



(eqn 3.56) 



i = V (i -P4 2 ). 



(eqn 3.57) 



The previously defined matrix with substituted values as calculated above satisfies the 
required property that P'P = 

D. EFFECTIVENESS OF THE AR(1,4) PROCEDURE 

As with any new procedure, it is important that its effectiveness relative to older 
procedures be determined. In this case the AR(1,4) procedure will be compared to the 
two-step procedure. Specifically, these procedures shall be compared in three areas: 

• Accuracy of estimate for pq, 

• Ability to detect the AR(1,4) disturbance 

• Accuracy of estimates for regression parameters 

Two simulation programs written in APL were utilized to generate data for the 
analysis required to make the comparison. These programs are displayed in the 
appendix. Program Stats is used to generate data for the comparison of the estimates 
of pq and the ability to detect the AR(1,4) disturbance. This program utilizes the APL 
programs Norrand and Unirand for the generation of normally distributed random 
numbers. Norrand and Unirand are generic APL random number generators on the 
IBM 370/3033 AP computer system which was used to perform the simulations. After 
the program generates regression residuals, the Durbin-Watson, Wallis and Schmidt 
statistics are calculated. The Durbin-Watson and Schmidt statistics are then compared 
to their respective bounds for the selected number of observations and independent 
variables. If the statistic detects the disturbance, a counter will be incremented. In 
addition to determining whether a detection is made or not, the statistics are used in 
previously discussed formulas to calculate values for pj and pq. 
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Program Regress is the second simulation used. Its specific purpose is to 
generate data which will allow the comparison of regression parameters from both the 
two-step and AR(1,4) procedures. This program as does the previous, utilizes the 
Norrand and Unirand routines to generate normally distributed random numbers. In 
the program, values for one independent variable and intercept term were generated for 
a predetermined number of observations, from a normal distribution with zero mean 
and a variance of .0625. After generating the error term from the AR(1,4) model for 
specified values of Pj and p^, the dependent variables are determined by way of the 
general linear model with a value for P of one. Using this data, estimates of P are 
determined by way of ordinary least squares, the two-step procedure and the AR(1,4) 
procedure. 

As previously discussed in Chapter I, it was decided to determine how well the 
AR(1,4) procedure performs relative to the two-step procedure for a best case and 
worst case scenario. For the best case scenario, the number of observations was 100 
and the error term was generated from a normal distribution with zero mean and a 
variance of .01. For the worst case the number of observations was 20 and the 
variance of the normal distribution was 10. For both the best and worst case 
scenarios, 200 replications were performed for 17 combinations of Pj and p^. The 
number of combinations which could be investigated was limited by the initial terms of 
the P matrix which contained square roots. The problem was that for certain 
combinations of pj and p^ the term under the square root became negative. The 17 
combinations of Pj and p^ which were investigated were 
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14 P| — .7 p 4 — .3 

15 Pj = .7 p 4 = *5 

16 pj = .9 p 4 = .1 

17 pj = .9 p 4 = .3 . 

For a given value of p j , any value of p 4 greater than the displayed values for p 4 
caused the problem with terms under square roots in the P matrix. In subsequent 
graphs used to display results, each combination of Pj and p 4 will be referred to by the 
number as assigned above. Thus combination one is given where Pj = .1 and p 4 = 
.1. The effectiveness of the procedure in handling negative autocorrelation was not 
investigated. 

Mean square error (MSE) was selected to be the measure by which the best 
procedure would be selected. In the comparison of^s from both procedures, MSE 
was estimated separately for the slope and intercept terms utilizing the formula 




(eqn 3.58) 



where u was the actual value of the slope or intercept and n was the number of 
simulation replications. Specifically, the value of u was unity for this application. In a 
similar manner, the \fsfe for the vector p consisting of the slope and intercept term 
was determined via the formula 

MSE = £ ((f 0 - l)(fj - l)) 2 /n (eqn 3.59) 

where ^ is the intercept term and Pj is the slope term. Mean square error for the 
comparison of p 4 's was calculated in a similar manner. 

Power curves were constructed to allow the determination of which procedure 
could best detect the AR(1,4) disturbance. These curves were constructed by plotting 
for each Pjand p 4 combination the number of times each procedure detected the 
disturbance out of 100 possible attempts. As previously mentioned, if a trial turned 
out to bt inconclusive, it was not counted as a detection. The number of times for 
each Pj.p 4 combination that the procedures were inconclusive was plotted to enable 
readers who consider inconclusive trials as detections to reevaluate the power curves. 
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E. RESULTS 

Comparisons of mean square error for estimates of are shown in Figures 3.1 
and 3.2. Figure 3.1 displays mean square errors for the worst case scenario while 
Figure 3.2 displays MSEs for the best case situation. Mean squared errors for both the 
best and worst case scenario are consistent for both the AR(1,4) and two-step 
procedure estimates. For the best case scenario the two-step estimate is clearly 
superior for Pj,p^ combinations greater than nine. MSEs for two-step and AR(1,4) 
procedure estimators for the worst case scenario diverge for p j ,p^ combinations greater 
than 15. This indicates that for the worst case scenario, the AR(1,4) procedure 
estimate for p^ can be expected to perform as well as the two-step estimate. 

As previously discussed, an individual using the two-step procedure to test for an 
AR(1,4) disturbance might initially use the Durbin-Watson statistic. If the test 
indicated that autocorrelation did not exist, further testing with the Wallis statistic 
might not be considered. Since this was the case assumed, power curves were created 
which compared the Schmidt and Durbin-Watson statistics. The Schmidt statistic is 
considered as part of the AR(1,4) procedure while the Durbin-Watson is part of the 
two-step procedure. For the best case scenario with T = 100 and x = .01, shown in 
Figure 3.4, the Schmidt statistic except for the sixth combination is more powerful than 
the Durbin-Watson statistic. For this scenario, both statistics are equally powerful for 
Pj.Pq combinations greater than 10. For the worst case situation, shown in Figure 3.3, 
both statistics are equally capable across all combinations in detecting the AR(1,4) 
disturbance. As previously discussed the power curves for both procedures do not 
include inconclusive test results. Figures 3.5 and 3.6 display the number of times that 
the test statistics were inconclusive for the worst and best case scenarios respectively. 

Figures 3.7 through 3.12 show comparisons of mean square errors for the 

regression parameters of the two-step and AR(1,4) procedures. Figures 3.7 and 3.8 

show \ISEs calculated using both Pq and Pj for worst and best case situations, 

respectively. For the worst case with T = 20, x = 10 scenario, shown in Figure 3.7, 
A 

MSEs for both procedures follow similar trends which are consistent with each other. 
Figures 3.9 and 3.10 which display MSEs for the intercept and slope coefficients, 
respectively, show trends which support the previous findings of Figure 3.7. Mean 
square errors for the best case scenario are shown in Figures 3.8, 3.11 and 3.12. As 
shown in Figure 3.8, both procedures follow similar trends for Pj.P/j combinations up 
to 15. Clearly the two-step procedure diverges at the Pj.Pq combination 15. Figures 
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Figure 3.1 MSB Comparison of Estimators for p_j Worst Case 




Figure 3.2 MSE Comparison of Estimators for Best Case 



45 




Figure 3.3 Power of Statistics Worst Case 




Figure 3.4 Power of Statistics Best Case 
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Figure 3.5 Inconclusiveness of Statistics Worst Case 




Figure 3.6 Inconclusiveness of Statistics Best Case 
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3.11 and 3.12 indicate that this divergence is related to a divergence of the intercept 
estimate for the two-step procedure. Overall, it appears that regression parameters 
from both procedures are consistent with each other and that one procedure is not 
superior to the other with regards to the estimation of regression parameters. 

F. CONCLUSIONS AND RECOMMENDATIONS 

Overall, the AR(1,4) procedure as previously developed is just as capable as the 
two-step procedure in correcting the AR(1,4) disturbance for the Pj.p^j combinations 
for which it can be used. The Schmidt statistic as part of this procedure is effective in 
detecting the AR(1,4) disturbance. In certain situations, its abilities surpass those of 
the Durbin-Watson statistic from the two-step procedure. In particular, its capabilities 
arc noteworthy at higher values of T, lower values of and lower p j ,p^ combinations. 
It appears from the results that estimates of Pj calculated via the two-step procedure 
are better than estimates from the AR(1,4) procedure, at least for the order of 
calculation of p_j which was tested. If an individual suspects that an AR(1,4) 
disturbance is present in the data then it is recommended that the Schmidt, 
Durbin-Watson, and Wallis statistics be calculated. The Schmidt statistic should be 
utilized to detect the disturbance while the Durbin-Watson and Wallis statistics are 
used to calculate values of pj and p^. If the values for Pj and p_j are acceptable, then 
they should be used in the AR(1,4) P matrix to correct for the disturbance. If the pj 
and p^j values fall out of the acceptable range for the AR(1,4) procedure then the 
two-step procedure must be used. 

G. AREAS OF FUTURE RESEARCH 

Three areas of future research are suggested as a result of this research. The first 
of these would be to find a P matrix for the AR(1,4) procedure which was not limited 
by values of pj and p^. A second area of research would be to derive a correction 
procedure for a general fourth order autoregressive process. This derivation would 
follow closely what was done above for the AR(1,4) process. 

A final area of research would be to investigate the maximum likelihood 
estimation of a general AR(4) or an AR(1,4) process. Beach and MacKinnon [Ref. 11] 
have investigated the second order process and suggest possible ways to extend this 
procedure to the general fourth order case. 
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Figure 3.7 Comparison of Two-step and AR(1,4) Procedure's MSEs Worst Case 




✓N 

Figure 3.8 Comparison of Two-step and AR(1,4) Procedure's MSEs Best Case 
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Figure 3.9 Comparison of Two-step and AR(1,4) Procedure's Intercept MSEs Worst Case 




A 

Figure 3.10 Comparison of Two-step and AR(1,4) Procedure's Slope MSEs Worst Case 
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Figure 3.11 Comparison of Two-step and AR(1,4) Procedure's Intercept MSEs Best Case 
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Figure 3.12 Comparison of Two-step and AR(1,4) Procedure's Slope MSEs Best Case 
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APPENDIX 
PROGRAM LISTINGS 



C 

c 



THIS PROGRAM COMPUTES THE LOWER BOUND FOR THE SCHMIDT STATISTIC 
USING IMHOFF'S INTEGRAL. THE FILE MAT IS A FILE OF EIGENVALUES 
OBTAINED FROM THE SUBROUTINE LATRV. THE USER IS REQUIRED TO SET 
THE NUMBER OF OBSERVATIONS AND THE NUMBER OF INDEPENDENT VARIABLES. 
ADDITIONALLY, A STARTING FOR XKRITL SHOULD BE SUPPLIED. 

THE PROGRAM WILL PRINT VALUES FOR THE INTERGRAL AND XKRITL 
AS THE PROGRAM EXECUTES. IN THIS WAY THE USER IS ABLE TO 
SEE THE HOW WELL THE INTEGRAL IS CONVERGING. IF THE INTEGRAL 
IS CONVERGING TOO SLOWLY THEN THE VALUE OF XKRITL SHOULD BE 
CHANGED 

^A A tt^? 0 ^^afFg%«} 0 L?®i^ 100) ' L(100) ’ U(100) ' S(100) 

INTEGER IND.N.iS.K.t.Nft.KP 

open(unit=i6 , Eil£= mAt 1 1 

OPEN(UNIT=12 FILE='OUT') 

XKRITL=0. 

X K R I T L= 3 . 40 

ASSIGN VALUE TO ORDER OF MATRIX 
N=100 

SET NUMBER OF INDEPENDENT VARIABLES (INCLUDE INTERCEPT TERM 
K- 6 

IS=N*N 

SET DESIRED TRUNCATION ERROR 
EPS 1=. 0001 

SET DESIRED ACCURACY OF NUMERICAL INTEGRATION 
EPS2=. 0001 

SET DESIRED ACCURACY OF DISTRIBUTION FUNCTION 
EPS3=. 001 

SET VALUE OF DISTRIBUTION FUNCTION 

FD=. 05 

SET INPUT CODE 

IND=1 

ADJUST FOR THE NUMBER OF COLUMNS OF THE MATRICES 'A' AND 
Y X‘ WHICH ARE LINEAR COMBINATIONS OF EACH OTHER 
T=N-K 

GENERATE THE 'A' MATRIX 



C 

C 



CALL SCHMDT 
WRITEflO 
REWIND lfl 
READ 



,J), J=1,N), 1=1, N) 



IN 



READ 



'A' MATRIX FOR WHICH EIGENVALUE AND EIGENVECTORS 
WILL BE COMPUTED 

omiWe ( tiSe Eigenvalues and eigenvectors of matrix 'a 1 



(10 * 

c 
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CALL LATRV 
WRITEC 12 , * 

REWIND 12 
READ( 12 *_)(D( I) . 1=1, N) 

SET THE NUMBER dF non-zero eigenvalues 



NR=T 

SORT THE 
KP=K-1 
DO 10 J=1 ,T 

W J+KP) 



N-K SMALLEST EIGENVALUES 



10 CONTINUE 

COMPUTE LOWER BOUND FOR STATISTIC 
CALL FQUADfL, NR. FD, EPS1.EPS2 , XKRITL) 
XKRITD=ABS(XKRltL-tfKRITd) 
XKRITO=XKRiTL 

IF((ABS(FD-.05)-EPS3). LE.O) GO TO 30 
IF((FD-.05). LT. O) THEN 
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30 



XKRITL=XKRITL+. 5*XKRITD 
PRINT * FD.XKRITL 
GO TO 46 
END IF 

IF((FD-.05).GT. 0) THEN 
XKRITL=XKRITL-. 5*XKRITD 
PRINT * FD,XKRITL 
GO TO 46 
END IF 

PRINT *,FD,XKRITL 
STOP 
END 

SUBROUTINE LATRV 



PURPOSE: 

MATRIX 



FIND THE EIGENVALUES AND VECTORS OF A REAL SYMMETRIC 



DESCRIPTION OF VARIABLES: 

R= A REAL SYMMETRIC MATRIX, DESTROYED DURING COMPUTATIONS 
AND REPLACED BY THE TRANSPOSED MATRIX OF EIGENVALUES 
N= ORDER OF THE MATRICES 
V= MATRIX OF EIGENVECTORS 
D= VECTOR OF EIGENVALUES IN DESCENDING ORDER 
IND= INPUT CODE 

IF IND=0 EIGENVALUES AND EIGENVECTORS ARE 
COMPUTED 

IF IND=1 EIGENVALUES ARE COMPUTED ONLY 

METHOD: DIAGONALIZATION METHOD OF JACOBI 

REFERENCE: NATIONAL PHYSICAL LABORATORIES, 'MODERN COMPUTING 
METHODS^ 



SUBROUTINE LATRV (R,N, V.D, IND) 
DIMENSION R( 1) ,V( 1 ) , D( 1 J 
INDEX(I,J)=I+/6-rrN 

GENERATE IDENTITY MATRIX 
IF(IND) 1,1,4 

1 NN=N*N 

IJ=N+1 
DO 2 1=1, NN 

2 V(I)=0 

DO 3 1=1 , NN , I J 

3 V(I)=1 

COMPUTE FINAL NORM EPS 

4 EPS=0 

DO 6 1=1, N 
DO 6 J=1 N 
IJ=INDEX(j ,1) 

IF(I-J) 5.6.5 

5 EPS=EPS+R( lj)**2 

6 CONTINUE 

I F( EPS-. 0000001) 95,95.7 

7 EPS=SQRT(EPS)*. 0000061/^ 

DREMP=1 

8 L2=0 

DO 80 LR=1,N 
IR=INDEX(0’LR) 

DO 80 LK=Lft,N 
IF(LK-LR) 86,80,10 
10 IK=INDEX(6,LK) 

K1=LR+IK 

K2=LR+IR 

K3=LK+IK 

IF(ABS(ARK)-DREMP) 80,80,20 
20 L2=l 

COMPUTE SIN AND COS 
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c 

c 




25 

26 



;lr,j) 

LK,J) 



40 



45 



50 

80 

90 



ARR=F 
AKK= 

XL= ' 

XM=ARR-AKK 
SQ=SQRT(XL*XL+XM*XM) 

S IN2T=XL/SQ 
IF(XM) 25.26,26 

$in2t=-sin£t , x w x 

/Sq)/2) 

SINT2=SINT**2 
COST2=COST**2 

MODIFY MATRICES 
DO 40 J=1,N 
K4=J+IR 
K5=J+IK 
K6=INDEX( 

K7=INDEX 
ARJ=R(K4 
AKJ=R(K5 

R( K4)-ARj*COST+AKJ*SINT 
R(K5)=AKJ*COST-ARJ*SINT 
R(K6)=R(K4) 

R(K7)=R(K5) 

CONTINUE 
K7=LK+IR 
R(K1)=0 
K7 )=0 

'K2)=ARR*COST2+AKK*SINT2+ARK*SIN2T 
' K3)=ARR*SINT2+AKK*COST2-ARK*SIN2T 
: (IND) 45,45,80 
__ 50 1=1, N 
K1=I+IR 
K2=I+IK 
ARR=V(K1) 

AKK=V(K2) 

V(K1)=C0$T*ARR+SINT*AKK 
V( K2)-COST*AKK-SINT*ARR 
CONTINUE 
IF(L21 8.90,8 
DREMP=DREfaP*. 1 

COMPARE THRESHOLD 



DREMP WITH NORM EPS 



95 

100 

110 

115 



120 

130 



135 

140 



IF(DREMP-EPS) 95,8,8 

SORT EIGENVALUES AND EIGENVECTORS 



DO 100 1=1. N 
I J= INDEX ( f , I ) 

D(I )=R( IJ; 

DO 130 1=1, N 
DO 130 J=I,N 
IF(D(J)-D(1)) 
ARR=D(I) 

D( I )=D( J ) 

D(j1=ARR 

IF(IND) 115,115,130 
DO 120 K=1,N 



130,130,110 



lJ=INDEX(fc,r 
ARK=v(rj; 



IK=INDEX(K;j: 



vru }=v(iK) 

V(IK)=ARK 

CONTINUE 

CONTINUE 

REPLACE R BY TRANSPOSED 
IF(IND) 135,135,150 
DO 140 1=1, N 
DO 140 J=1,N 
IJ=INDEX(J I) 

IK=INDEX(I J) 

R( I J )=V( IK) ' 



MATRIX OF EIGENVECTORS 
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150 



RETURN 

END 



SUBROUTINE FQUAD 



. LE. X) OF 



PURPOSE: FIND THE DISTRIBUTION FUNCTION F(X)=P(Q 
QUADRATIC FORMS IN NORMAL VARIABLES. 

DESCRIPTION OF PARAMETERS 

S = A VECTOR OF LENGTH NR CONTAINING THE NON-ZERO 
EIGENVALUES 

NR = NUMBER OF NON-ZERO EIGENVALUES 

FD = VALUE OF THE DISTRIBUTION FUNCTION 

EPS1 = DESIRED TRUNCATION ERROR TU 

EPS2 = DESIRED ACCURACY OF NUMERICAL INTEGRATION 

XKRIT = VALUE OF X 

METHOD: THE INTEGRAL DERIVED BY IMHOF IS USED AND NUMERICALLY 
INTEGRATED BY SIMPSON T S RULE 

REFERENCE: 'COMPUTING THE DISTRIBUTION OF QUADRATIC FORMS IN 
NORMAL VARIABLES' BY J. P. IMHOF 
BIOMETRIKA (1961), 48, 3 AND 4, P.419 



C 

c 



F^UAD(S 



SUBROUTINE 
DIMENSION S( 
FD=0. 

UB=0. 

C2=0. 

C 1=0. 

FBL=0. 

VINT2=0. 

FD1=0. 

SUMK=0. 

XK=0. 

Si_AM=0. 

DO 6 1=1, NR 
S^=S(.)-X KRIT 



, NR, FD, EPS1 , EPS2 , XKRIT) 



SLAM COMPUTE 5 U A P L P^4 B 0 S ^W ^INTEGRAL GIVEN THE TRUNCATION 



ERROR. 



UB=EXP(-. AL0G(EPS1*XK) + 1. 
PRINT ,UB / . 



TU=EXP( 1. 14472989 
PRINT *,TU 
VA.=1.>TU-EPS1 
PRINT * , VAL 
IF((1. /'TU-EPSl) 

I FCC TU— EPS 1 ) . EQ 
~~ l./TU-EPf' 



14472989 + SLAM)/XK) 
ALOG(XK) + SLAM + XK*ALOG(UB)) 



LT 



SI). G 






0) GO TO 20 
GO TO 9 
0) GO TO 9 



10 

11 



14 



15 

16 



IF( H . 

UB=UB + 5. /XK 
GO TO 7 

COMPUTATION OF THE INTEGRAND GIVEN THE VALUE OF U 
I F(U) 15,15,11 
TETA=Q. 

RH0=1. 

DO 14 1=1, NR 
C1=S( I )*U 

C2=l. + Cl**2 / x 

T ETA=TETA + ATAN(Cl) 

RH0=RH0*(C2**. 25) 

T ETA=. 5^(TETA) 

FBL=SIN(TETA)7 



GO TO 18 
FBL=0. 

DO 16 1=1 



ETA)/(RHO*U) 



FBL=F3L+SU) 
FBL=. 5*( FBL 



NR 



■XKRIT) 
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18 GO TO (21.22,23.24,25). KSKIP 

EVALUAYkJn OF fHE TNTEGRAL BY SIMPSON'S RULE 

20 RANGE=UB 

MH=1 

U=RANGE*. 5 
KSKIP=1 
GO TO 10 

21 SUMK=FBL*RANGE*2./3. 

U=0. 

KSKIP=2 
GO TO 10 

22 V I NT2=SUMK+FBL* RANGE/6. 

U=UB 
KSKIP=3 
GO TO 10 

23 VINT2=VINT2 + FBL*RANGE/6. 

FD=. 5-. 318309886*VINT2 
DO 28 NIT=1 , 14 
FD1=FD 

VINT2=(VINT2-SUMK*. 5)*. 5 

MH=2*MH 

STEP=RANGE/MH 

U=STEP*. 5 

KSKIP=4 

GO TO 10 

24 SUMK=FBL 

IF(MH. EQ. 2) THEN 

U=U+STEP 

KSKIP=5 

GO TO 10 

END IF 

DO 25 K=2 ,MH 
U=U+STEP 
KSKIP=5 
GO TO 10 
SUMK=SUMK+FBL 
SUMK=SUMK*STEP*2./3. 

V I NT 2=V I NT 2 + SUMK 
FD=. 5 - ,318309886*VINT2 
IF(NIT-3) 28,28.27 
IF(ABS(FDl-FD)-Ef > S2) 29,28,28 
CONTINUE 
PAUSE 1234 
RETURN 
END 

SUBROUTINE SCHMDT 
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28 
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PURPOSE: GENERATE THE 
COMPUTED 



'A' MATRIX FOR WHICH EIGENVALUES ARE 



DESCRIPTION OF PARAMETERS 

A = MATRIX FOR WHICH EIGENVALUES ARE COMPUTED 
N = ORDER OF MATRIX 



SUBROUTINE SCHMDT(A,N) 
REAL A( 100 , 100) 

DO 10 1=1, N 
DO 20 J=1 N 
A(I.J)=0. 

20 CONTINUE 
10 CONTINUE 
A 1,1)=2. 

Af 1 2)=-l. 

A( 1 5 =-l. 

DO 3o 1=2 , N 
J=I 

IFfI.LE.4) THEN 
A(I,J)=3. 
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;j+n)=-i. 



IF(( I. GE. (N-3) ). AND. (I. LT. N)) THEN 

A( I ,J )=3. 

a( i ; ( j+i) )=-i. 



30 




IF(I.EQ.N) THEN 
I )-2, 



4). AND. (I. LT. (N-3))) THEN 



CONTINUE 

RETURN 

END 
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THIS PROGRAM COMPUTES THE UPPER BOUND FOR THE SCHMIDT STATISTIC 



LE MAT IS A FILE OF EIGENVALUES 
RV. THE USER IS REQUIRED TO SET 
HE NUMBER OF INDEPENDENT VARIABLES. 



C 

C 



C 

C 



OPEN( _ _ 
OPENfUNIl 
XKRIT ‘ 



USING IMHOFF'S INTEGRAL THE F 
OBTAINED FROM THE SUBROUTINE LA 

THE NUMBER OF OBSERVATIONS AND .... 

ADDITIONALLY, A STARTING FOR XKRITU SHOULD BE SUPPLIED. 

THE PROGRAM WILL PRINT VALUES FOR THE INTERGRAL AND XKRITU 
AS THE PROGRAM EXECUTES. IN THIS WAY THE USER IS ABLE TO 
SEE THE HOW WELL THE INTEGRAL IS CONVERGING. IF THE INTEGRAL 
IS CONVERGING TOO SLOWLY THEN THE VALUE OF XKRITU SHOULD BE 
CHANGED 

m M°MHr o) ' D(ioo) ' L(loo) - u(ioo) ' s(ioo) 

INTEGER IND.N.lS.K.t.Nft.KP 
.unit=iMil£='mAt') 

UNIT=12 FILE='OUT') 

ru=o. 

XKRITU=3 80 

ASSiGN VALUE TO ORDER OF MATRIX 

N=20 

SET NUMBER OF INDEPENDENT VARIABLES (INCLUDE INTERCEPT TERM 
K=5 

j ^ ^ j^j 

SET DESIRED TRUNCATION ERROR 
EPS1=. 0001 

SET DESIRED ACCURACY OF NUMERICAL INTEGRATION 
EPS2= 0001 

SET DESIRED ACCURACY OF DISTRIBUTION FUNCTION 
EPS3=. 001 

SET VALUE OF DISTRIBUTION FUNCTION 

FD=. 05 

SET INPUT CODE 

IND=1 

ADJUST FOR THE NUMBER OF COLUMNS OF THE MATRICES 'A' AND 
Y X' WHICH ARE LINEAR COMBINATIONS OF EACH OTHER 
T=N-K 

GENERATE THE 'A' MATRIX 



CALL SCHMDT 
WRITE! 10 *' 
REWIND 1(5 
READ IN 



,J), J=1,N), 1=1, N) 



■A 1 MATRIX FOR WHICH EIGENVALUE AND EIGENVECTORS 
WILL BE COMPUTED 
READ (10 *5 (R(I). 1=1, IS) 

OMI^UTE THE EIGENVALUES AND EIGENVECTORS OF MATRIX 'A' 
CALL LATRVfR.N.V.D.IND) 



40 REWIND 12 

READ( 12 *_)(D( I) . 1=1. N) 

C SET THE NUMBER dF NON-ZERO EIGENVALUES 

NR=T 

C SORT THE N-K LARGEST EIGENVALUES 

DO 20 1=1, T 

um=Dm 

20 CONTINUE 

C COMPUTE LOWER BOUND FOR STATISTIC 

CALL FQUADTu. NR, FD.EPS1.EPS2. XKRITU) 

C SUCCESSIVELY HAEf XKRITU fO FIND THE 5 PERCENT 

C , SIGNIFICANCE POINT 

XKRITD=ABS(XKRITU-XKRITO) 

XKRITO=XKRITU 

IF(!ABS(FD-. 05)-EPS3). LE. 0) GO TO 30 
IF((FD-. 05). LT. O) THEN 
XKRITU=XKRITU+. 5*XKRITD 
PRINT *.FD, XKRITU 
GO TO 46 



END IF 
IF((FD 



xkpiiu=XkrKu t ‘ q) 



PRINT 
GO TO 



„ THEN 
5*XKRITD 
*.FD, XKRITU 
46 
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END IF 

30 PRINT *,FD,XKRITU 
STOP 
END 

SUBROUTINE LATRV 



PURPOSE: 

MATRIX 



FIND THE EIGENVALUES AND VECTORS OF A REAL SYMMETRIC 



DESCRIPTION OF VARIABLES: 

R= A REAL SYMMETRIC MATRIX, DESTROYED DURING COMPUTATIONS 
AND REPLACED BY THE TRANSPOSED MATRIX OF EIGENVALUES 
N= ORDER OF THE MATRICES 
V= MATRIX OF EIGENVECTORS 
D= VECTOR OF EIGENVALUES IN DESCENDING ORDER 
IND= INPUT CODE 

IF IND=0 EIGENVALUES AND EIGENVECTORS ARE 
COMPUTED 

IF IND=1 EIGENVALUES ARE COMPUTED ONLY 

METHOD: DIAGONALIZATION METHOD OF JACOBI 

REFERENCE: NATIONAL PHYSICAL LABORATORIES, 'MODERN COMPUTING 
METHODS' 



SUBROUTINE LATRV (R,N,V.D,IND) 
DIMENSION R(1 ) . V( 1) ,D( l) 
INDEXa.J)=I+fJ-lJ*N 

GENERATE IDENTITY MATRIX 
IF(IND) 1,1,4 

1 NN=N*N 

IJ=N+1 
DO 2 1=1, NN 

2 V(I)=0 

DO 3 1=1 ,NN , IJ 

3 V(I)=1 

COMPUTE FINAL NORM EPS 

4 EPS=0 

DO 6 1=1, N 
DO 6 J=1,N 
IJ=INDEX(J,I) 

IF(I-J) 5,6,5 

5 EPS=EPS+R(lj)**2 

6 CONTINUE 

IF( EPS-. 0000001) 95,95.7 

7 EPS-SQRT(EPS)*. 00000t)l/N 

DREMP=1 

8 L2=0 

DO 80 LR=1,N 

ir=index(o!lr) 

DO 80 LK=LA a N 
IF(LK-LR) 86,80,10 
10 IK=INDEX(6,LK) 

K1=LR+IK 

K2=LR+IR 

K3=LK+IK 

IF(ABS(ARK)-DREMP) 80,80,20 
20 L2=l 

COMPUTE SIN AND COS 
ARR=R(K2) 

AKK=R(K3) 

XL=2*ARK 
XM=ARR-AKK 
SQ=SQRT(XL*XL+XM*XM) 
SIN2T=xL/SQ 
IFCXM) 25,26,26 
25 SIN2T=-SIN^T 
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C 

C 



;lr,j) 

LK,J) 



40 



45 



50 

80 

90 



95 

100 

110 

115 



120 

130 



135 



140 

150 



« r n t ^^s x W /SQ)/2) 

SINT2=SINT**2 
COST2=COST**2 

MODIFY MATRICES 
DO 40 J=1,N 
K4=J+IR 
K5=J+IK 
K6=INDEX( 

K7= I NDEX ( 

ARJ=R(K4 
AKJ=R(K5! 

R( K4)-ARJ*C0ST+AKJ*SINT 
R(K5)=AKJ*C0ST-ARJ*SINT 
R(K6)=R(K4) 

R(K7)=R(K5) 

CONTINUE 

K7=LK+IR 

R(K1)=0 

R(K7}=0 

R( K2)=ARR*C0ST2+AKK*SINT2+ARK*SIN2T 
RLK3)=ARR*SINT2+AKK*COST2-ARK*SIN2T 
IF(IND) 45,45,80 
DO 50 1=1, N 
K1=I+IR 
K2=I+IK 
ARR=V(K1) 

AKK=V(K2 1 

V(K1)=C0ST*ARR+SINT*AKK 
V(K2l=C0ST*AKK-SINT*ARR 
CONTINUE 
IF(L21 8.90,8 

dremp=6reMp*: 1 

COMPARE THRESHOLD DREMP WITH NORM EPS 
IF(DREMP-EPS) 95.8,8 

SORT EIGENVALUES AND EIGENVECTORS 
DO 100 1=1. N 
IJ=INDEX(f,I) 

D(I)=R(IJ) 

DO 130 1=1, N 
DO 130 J=I,N 

IF(D(J)-D(i)) 130,130,110 
ARR=D( I) 

D( I )-D( J) 

Dfj1=ARR 

IF(IND) 115,115,130 
DO 120 K=1 ,N 
IJ=INDEX(K,I) 

IK=INDEX(K J) 

ARK=V( IJ) 

V( IJ)=V( IK) 

v(ik1=ark 

CONTINUE 
CONTINUE 

REPLACE R BY TRANSPOSED 
IF(IND) 135,135,150 
DO 140 1=1. N 
DO 140 J=1,N 
IJ=INDEX(J I) 

IK=INDEX(I J) 

R(IJ)=VCIK) 

RETURN 
END 

SUBROUTINE FQUAD 



MATRIX OF EIGENVECTORS 



PURPOSE: FIND THE DISTRIBUTION FUNCTION F(X)=P(Q . LE. X) OF 
QUADRATIC FORMS IN NORMAL VARIABLES. 

DESCRIPTION OF PARAMETERS 
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S = A VECTOR OF LENGTH NR CONTAINING THE NON-ZERO 
EIGENVALUES 

NR = NUMBER OF NON-ZERO EIGENVALUES 

FD = VALUE OF THE DISTRIBUTION FUNCTION 

EPS1 = DESIRED TRUNCATION ERROR TU 

EPS2 = DESIRED ACCURACY OF NUMERICAL INTEGRATION 

XKRIT = VALUE OF X 

METHOD: THE INTEGRAL DERIVED BY IMHOF IS USED AND NUMERICALLY 
INTEGRATED BY SIMPSON'S RULE 

REFERENCE: 'COMPUTING THE DISTRIBUTION OF QUADRATIC FORMS IN 
NORMAL VARIABLES' BY J. P. IMHOF 
BIOMETRIKA (1961), 48, 3 AND 4, P.419 



C 

C 



SUBROUTINE FQUAD(S , NR , FD, EPS1 , EPS2 , XKRIT) 

DIMENSION S(l) 

FD=0. 

UB=0. 

C2=0. 

C1=0. 

FBL=0. 

VINT 2=0 . 

FD1=0. 

SUMK=0. 

XK=0. 

SLAM=0. 

DO 6 1=1, NR 

XK=Ik+ ( 5 _XKRIT 

SLAM=SLAM+. 5*ALOG(ABS(S( I))) 

COMPUTE UPPER-BOUND OF INTEGRAL GIVEN THE TRUNCATION 
ERROR. 

UB=EXP(-(AL0G(EPS1*XK) + 1.14472989 + SLAM)/XK) 

PRINT * ,Ut 



10 

11 



14 



15 

16 
18 
20 



21 



TU=EXP( 1. 14472989 
PRINT *,TU 
VAL=1. /TU-EPS1 
PRINT *,VAL 
IF( ( 1. /TU-EPS1). LT. 0) 
IFC(TU-EPSl). EQ. 0} GO 
I Ff Cl- /TU-EPS1). GT. 0) 
UB=UB + 5./XK 
GO TO 7 

COMPUTATION OF THE 
I F(U) 15,15,11 
TEl A=0. 

RHO=l. 

DO 14 1=1, NR 
C1=S(I)*U 

C2=l. + C 1**2 , v 

TETA=TETA + ATAN(Cl) 
RHO=RHO*(C2**. 25) 

TETA=. 5*(TETA) 

FBL=SIN( TETA)/ 



ALOG(XK) + SLAM + XK*ALOG(UB)) 



GO TO 
TO 9 
GO TO 



20 



INTEGRAND GIVEN THE VALUE OF U 



18 



ETA)/(RHO*U) 



GO TO 
FBL=0. 

DO 16 1=1. NR 
FBL=FBL+SU) 

FBL=. 5*f FBL— XKRIT ) 

GO TO (21.22,23,24:25). KSKIP 

EVALUXTlQN CfF fHE INTEGRAL BY SIMPSON'S 
RANGE=UB 
MH=1 

U=RANGE*. 5 
KSKI P=1 
GO TO 10 

SUMK=FBL*RANGE*2. /3. 



RULE 
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u=o. 

KSKIP=2 
GO TO 10 

22 VINT2=SUMK+FBL*RANGE/6. 

U=UB 
KSKIP=3 
GO TO 10 

23 VINT2=VINT2 + FBL*RANGE/6. 

FD=. 5-. 318309886*VINT2 
DO 28 N I T= 1 , 14 
FD1=FD 

VINT2=(VINT2-SUMK*. 5)*. 5 

MH=2*MH 

STEP=RANGE/MH 

U=STEP*. 5 

KSKIP=4 

GO TO 10 

24 SUMK=FBL 

IFCMH. EQ.2) THEN 

U=U+STEP 

KSKIP=5 

GO TO 10 

END IF 

DO 25 K=2 , MH 
U=U+STEP 
KSKIP=5 
GO TO 10 
SUMK=SUMK+FBL 
SUMK=SUMK*STEP*2. /3. 
VINT2=VINT2 + SUMK 
FD=. 5 - , 318309886*VINT2 
I FC NIT— 3^ 28,28,27 
IF(ABS(FD1-FD)-ePs2) 29,28,28 
CONTINUE 
PAUSE 1234 
RETURN 
END 

SUBROUTINE SCHMDT 
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PURPOSE: GENERATE THE 
COMPUTED 



'A' MATRIX FOR WHICH EIGENVALUES ARE 



DESCRIPTION OF PARAMETERS 

A = MATRIX FOR WHICH EIGENVALUES ARE COMPUTED 
N = ORDER OF MATRIX 



SUBROUTINE SCHMDT(A,N) 
REAL A(IOO.IOO) 

DO 10 1=1, N 
DO 20 J=1 ,N 
A(I.J)=0. 

20 CONTINUE 
10 CONTINUE 
A( 1,1 )=2. 

Afl 2)=-l. 

Afl’5)=-1. 

DO 50 I=2,N 

IFfl.LE. 4) THEN 

Au,gS=3 ' 

=-i: 

=-i. 




AND. (I. LT.N)) THEN 
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Stf-K- 4 ”- 1 - 

IF(I.EQ.N) THEN 
IF(( I. GT. 4). AND. (I. LT. (N-3))) THEN 




30 



RETURN 
END 
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Program Stats is an APL simulation which was used to generate data for the 
comparison of the estimates of and the ability of the Schmidt statistic as part of the 
AR(1,4) procedure and the Durbin- Watson statistic as part of the two-step procedure 
to detect the AR(1,4) disturbance. 

V STATS ;ADliAD;AWlBi AW1 1 AW2 ; AW 3 ; AWnB ; AW*+ ; AW5 ; AW ; AS ; PW 
[ 1 ] n ENTER VALUES FOR RHO 1 AND RHOn 
[ 2 ] ' ENTER A VALUE FOR RHO 1 ' 

[3] RHOl+O 

[ 4 ] ' ENTER A VALUE FOR RHO 4 ' 

[5] RHOV+Q 

[ 6 ] ' ENTER NUMBER OF OBSERVATIONS • 

[7] NOBS * □ 

[ 8 ] ' ENTER VARIANCE OF NORMAL DISTRIBUTION ' 

[9] SIG2+U 

[10] A INITIALIZE BOUND FOR STATISTICAL TESTS 

[11] DWL+l. 65 

[12] DWU+l . 69 

[13] DWAL*-1 . 59 

[14] DWAU+1 . 63 

[15] DSMTL+3 .45 

[16] DSMTU+3. 53 

[ 1 7] fi INITIALIZE OUTER LOOP VECTORS AND COUNTERS 

[18] RHOAV l-e-40p0 

[19] i?ff0AP4S«-4OpO 

[20] RHOAV^+HOpO 

[21] Z?l«-40p 0 

[22] i?4S<-40 p 0 

[23] £4-e-40p0 

[24] DWCNT<rO 

[25] DSTCNT<rO 

[26] DWACNT+ 0 

[27] DWINCT^O 
[2 8] STINCT<rO 

[29] DWAICT+Q 

[30] n GENERATE DURBIN -WATSON 'A' MATRIX 
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[31] NO+NOBS 

[32] ADlB*-( ( (NO -2 )p0),”l,2,~l)),( {NO -2 )pO ) , ”l , 1 

[33] AD1<-1,~1, (( (NO* (.NO-2) ) + (NO-2) )pADlB 

[34] AD+(NO ,NO)pADl 

[3 5] o GENERATE THE WALLIS 'A' MATRIX 
[36] AWlB*-( ( (NO -5 )-3 )pO ) 

[3 7] AW 1C*- ( ( ( (NO -5 ) + l)p0),l,0,0,0,~l)) ,AW1B ) 

[3 8] AJ71-«-l ,0,0,0, - l,((((3 *N0 )+3 ) pAWlC 

[3 9] AW2*-( ( (NO -8 )*N0 ) — 9 ) p ( ( (NO -8 )p0), _ l,0,0,0,2,0,0,0,~l) 

[40] AW3*-~ 1 , 0,0, 0,2,0, 0,0, ~1 ,AW2 

[41] AWUB+ (((3xNO ) + 3 ) )p ((.((NO- 5)+l)p0), "1,0, 0,0,1) 

[42] AWU*-( ( (NO- 5 )-3 )p0 ) , ~1 , 0 , 0 , 0 , 1 ,AWVB 

[43] AW5*-AW1 ,AW3 , AJ74 

[44] AW*-(NO ,NO)pAW5 

[45] fl GENERATE THE SCHMIDT 'A' MATRIX 

[46] AS+AD+AW 

[47] A OUTE^LOOP FOR REPLICATION 

[48] COUNT 0+0 

[49] OUTER :COUNTO+COUNTO+ 1 

[5 0] fl INITIALIZE INNER LOOP VECTORS 

[51] RHOHAT1*-10 p 0 

[52] RHOHATnS*- lOpO 

[53] RHOHATn+lOpO 

[ 54] fl INNER LOOP TO CALCULATE RHO 1 , RHOn AND DETERMINE DETECTION 

[55] COUNT I +0 

[56] INNER: COUNTI+COUNTI+l 

[57] DETDW+0 

[58] DETSMT+o 

[59] DETWA+0 

[60] INSTCT*-0 

[61] INDUCT *- 0 

[62] INDWAT+0 

[63] fl GENERATE THE RANDOM ERROR 

[64] V*-NOBS NORRAND (0 , SI G2) 

[65] E+NOBSpO 
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[ 66 ] N+ 4 

[67] A CALCULATE CORRECTION TERMS 

[ 68 ] e7-<-( 1 - (BB04* 2 ) )*0 . 5 

[69] I<r-(RH01*J) 

[70] G+RHOUxI 

[71] F+(1+(RH01*2 )- (RHOn*2 )- (J* 2 ) )*0 . 5 

[72] El<-(i?ff01*F) 

[73] D<— (I*G)*F 

[74] C+(1 + (.RH01*2 )- (RHOU*2 )- (Fl *2 ) )* 0 . 5 

[75] Bl-s-(Ftf01-(£xFl))*<; 

[76] A<r(l-(RHOU*2 )- (C *2 )-(£>*2 )-(Bl *2 ) )* 0 . 5 
[7 7] A CALCULATE THE FIRST FOUR ERROR TERMS 

[78] F[l]«-7[1]*4 

[79] E[2]<-(7[2]-(fllxff[l])) + C 

[80] F[ 3 ]«-( 7 [ 3 ]-(£xF[l] )-(FlxF[ 2 ] ))*F 

[81] F[4]«-(7[4]-(CxF[l] )-(JxF[3] ))*</ 

[82] MM:W«-W +1 

[83] F[W]«-(FB 01 xF[W-l] ) + (.RHOUxElN-Ul )+VlNl 

[84] +(N<NOBS)/MM 

[85] E+(NOBS,l)pE 

[ 86 ] A GENERATE THE INDEPENDENT VARIABLES 

[87] Xl+(NOBS , 1 )p (.NOBSpl ) 

[ 88 ] X2+(N0BS ,l)p (NOBS NORRAND 0 0.0625 ) 

[89] X«-X 1 ,X 2 

[ 9 0 ] A GENERATE THE TRUE BETA 
[91] BETAT+ 2 1 p 1 1 

[9 2 ] b GENERATE THE DEPENDENT VARIABLES 

[93] Y+(X+.xBETAT)+E 

[94] a LEAST SQUARES ESTIMATE OF BETA 

[9 5] B-e-YlX 

[96] A GENERATE THE RESIDUALS 

[97] mW-(X+.xfl) 

[98] A CALCULATE THE DURBIN -WATSON STATISTIC 

[99] DW<- , ( (§EHAT) + . x (AD+. *EHAT) )* ( (§EHAT) + . *EHAT) 
[10 0] A CALCULATE THE SCHMIDT STATISTIC 
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[101] DS+ , ( ( §EHAT ) + . X QS+ . *EHAT ) ) * ( ( §EHAT ) + . x£ff Ar ) 

[102] n COMPARE TABLED VALUE FOR DURBIN -WATSON STATISTIC TO 
[10 3] fi THE ABOVE CALCULATED VALUE FOR THE DURBIN -WAT SON 

[104] DETDW+iDWZDWL) 

[105] DWCNT+DWCNT+DETDW 

[106] INDWCT+ ( (DWL<DW )a (DWU>DW ) ) 

[107] DWINCT+DWINCT+INDWCT 

[10 8] fi COMPARE TABLED VALUE FOR SCHMIDT STATISTIC 

[109] A TO ABOVE CALCULATED VALUE FOR SCHMIDT STATISTIC 

[110] DETSMT+ ( DSZDSMTL ) 

[111] DSTCNT+DSTCNT+DETSMT 

[112] INSTCT+((DSMTL<DS')k(.DSMTU>DS)) 

[ 113] STINCT+STINCT+INSTCT 

[114] R CALCULATE THE DURBIN -WATSON ESTIMATE FOR RHOl 

[115] RHOHAT11COUNTI1+1- (DW*2) 

[116] fi CALCULATE THE SCHMIDT ESTIMATE FOR RHO^ 

[117] RHVS1+2-RHOHAT11COUNTI1 

[118] RHnS2+DS*2 

[119] RHOHAT^S [COUNT I] +RHHS1 -RHnS2 

[120] fi GENERATE TRANSFORM MATRIX FOR FIRST ORDER PROCESS 

[121] PD+, (l-(RHOHATlLCOUNTn*2)')*0.5 

[122] PW+( (NO*NO)-l )p ( ( (NO- 1 )p0 ) , (-RHOHAT1 [ COUNTI1 ) ,1 ) 

[123] PDW+(NO ,NO)pPD ,PW 

[124] fi TRANSFORM RESIDUALS 

[125] ESTAR+PDW+.xEHAT 

[126] fi CALCULATE THE WALLIS STATISTIC 

[127] DWA+ , ( ( §ESTAR ) + . x ( AW+ . *ESTAR ) ) * ( ( SESTAR ) + . xESTAR ) 
[12 8] fl COMPARE TABLED VALUE OF WALLIS STATISTIC 

[129] fl TO CALCULATED VALUE OF WALLIS STATISTIC 

[130] DETWA+(DWAZDWAL) 

[131] DWACNT+DWACNT+DETWA 

[13 2] INDWAT<r ( (DWAL<DWA)a(DWAU>DWA)) 

[133] DWAICT+DWAICT+INDWAT 

[134] a CALCULATE THE WALLIS ESTIMATE FOR RHOU 

[135] RHOHATmcOUNTIl+1- (DWA*2) 
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[136] +(C0UNTI<5 ) /INNER 
[13 7] RH0AV1 ZC0UNT01 + < +/RE0EAT1 )*5 
[13 8] RHOA V 45 [ CO UN TO ] <- ( + /RHORA T4S ) * 5 
[13 9] RE OA 74 ICO UNTO ]-«-(+ /RE OH A T4 ) t 5 

[140] +(COUNTO<nO ) /OUTER 

[141] CNTDW+DWCNT 

[142] CNTDST+DSTCNT 

[143] CNTDWA +DWA CNT 

[144] R1+RHOAV1 
[14 5] RHS+RHOAV^S 
[146] RU+RHOAVn 

V 
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Program Regress is an APL simulation which generates data for the comparison 
of the regression parameters from both the two-step and AR(1,4) procedures. 

V REGRESS -,AD iADUAWIB lAWl iAW2 ;AW3 iAWUB lAWn ;AW5 iAW iAS 
[ 1 ] A ENTER VALUES FOR RHO 1 AND RHOn 
[2 ] ' ENTER A VALUE FOR RHOl ' 

C 3 ] RHOl+U 

[4] ' ENTER A VALUE FOR RHOU ' 

[5] RHOU+U 

[ 6 ] ' ENTER NUMBER OF OBSERVATIONS ' 

[7] NOBS+U 

[ 8 ] 1 ENTER VARIANCE OF NORMAL DISTRIBUTION ' 

[9] SIG2+U 

[10] DSMTL+3.H 5 

[11] DSMTU+ 3.53 

[12] A GENERATE DURBIN -WATSON ' A ' MATRIX 

[13] NO+NOBS 

[14] ADlB+( ( (NO -2 )p0),~l,2,~l)), ( (NO-2 )p 0 ) , "l , 1 

[15] AD1+1,~1, (((NO*(NO-2))+(NO-2))pADlB 

[16] AD+(NO ,NO)pADl 

[17] A GENERATE THE WALLIS ' A ' MATRIX 

[18] AWlB+(((NO-5)-3)pO) 

[19] AWlC+( ( ( (.NO -5 )+l)p0),l,0,0,0,”l)) ,AW1B ) 

[20] ^J71-«-l ,0,0,0,”l,((( (3*N0 )+3 )pAWlC 

[21] AW 2+ ( ( (NO -3 )*NO ) — 9 ) p ( ( (NO -8 )p0),~l,0,0,0,2,0,0,0,~l) 

[22] AW3+~1, 0,0, 0,2, 0,0,0, ~1,AW2 

[23] AWUB+( ((3*N0 )+3 ) )p ( ( ( (W0-5 )+l )p0 ) , "l , 0,0, 0,1) 

[24] AWn<r(((NO-5)-3 )p0 >,"1,0,0,0,1,711/45 

[25] AW5+AW1 ,AW3 ,AWU 

[26] AW+(NO,NO)pAW5 

[27] A GENERATE THE SCHMIDT 'A' MATRIX 

[28] AS+AD+AW 

[ 2 9 ] A GENERATE TRANSFORM MATRIX FOR FIRST ORDER PROCESS 

[30] NO+NOBS 

[31] PD+,(1~(RH01*2))*0 .5 

[32] PW+((NO*NO)-l)p(((NO-l)pO) ,(~RHOl) ,1) 
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[33] PDW+(NO ,NO)pPD ,PW 

[34] n GENERATE TRANSFORM FOR FOURTH ORDER PROCESS 

[35] Pl«-(l-(tf//04*2))*0.5 

[36] P41«-P1, ((/VO-l)pO) 

[37] P42*-0,P1 , ( (NO-2 )pO ) 

[38] P43«-0 , 0 ,P1 , ( (NO -3 )p0 ) 

[39] P44«-0,0,0,P1, ((/V0"4)p0) 

[40] P4 5«- ( (/VO-4 )x/V0)p((-BB04), 0,0, 0,1, ((/VO-4 )p0)) 

[41] P4«-(JV0 ,/VO ) p (P41 ,P42 ,P4 3 ,P44 ,P4 5 ) 

[42] n INITIALIZE MATRICES AND VECTORS FOR OUTER LOOP 

[43] BLSIAV+VOpO 

[44] BLSSAV+^OpO 

[45] B2SIAV+U0p0 

[46] B2SS71V«-40p0 

[47] B14J717«-40p0 

[48] B14S717«-40p0 

[49] DSBCNT+0 

[50] DSACNT+ 0 

[51] STIBCT+ 0 

[52] STANCT+ 0 

[ 53 ] A OUTER LOOP FOR REPLICATION 

[54] COUNTO <- 0 

[55] OUTER : COUNTO+COUNTO + 1 

[ 56] p INITIALIZE MATRICES AND VECTORS FOR INNER LOOP 

[57] BLSI+ 5p0 

[58] BLSS+ 5p0 

[59] B2Si>5pO 

[60] B2SS+5pO 

[61] B14J«-5p0 

[62] P140^-5p 0 

[63] P INNER LOOP TO CALCULATE RHO 1 , RHOn AND DETERMINE DETECTION 

[64] COUNT 1 + 0 

[65] INNER : COUNTI+COUNTI+1 

[66] P GENERATE THE RANDOM ERROR 

[67] V+NOBS NORRANDiO ,SIG2) 
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[68] E+NOBS pO 

[69] N+ 4 

[70] o CALCULATE CORRECTION TERMS 

[71] - (RHOU*2 ) )*0 . 5 

[72] I-e--(Ptf01*c7) 

[73] C+RHOUxI 

[74] F«-(l + (m?l*2 )- (FF04*2 )- (1*2 ) )*0 . 5 

[75] E1+-(RH01*F') 

[76] D+-U*G)*F 

[7 7] C-e-(l + (FF01*2)-(FP04*2)-(Fl*2))*O.5 

[78] Bl*-(RH01-(DxEl))±C 

[79] 71-Kl - (FF04*2 )-(G*2 )- (D* 2 )-(B 1*2 ) )*0 . 5 

[ 8 0 ] R DETERMINE TRANSFORM MATRIX FOR 1 , 4 PROCESS 

[81] P141-e-d, ((iVO-l)pO) 

[82] P14 2<-B1 ,C,( (NO -2 )p0 ) 

[83] P143*P,F1,F, ( (iVO-3 )p0 ) 

[84] P144-«-(7, 0 ,1 ,J , ( (iVO-4 )p0 ) 

[8 5] P14 5<-((W-4)x^O)p((-PP04),0,0, ( -FP01 ) , 1 , ( (iV0-4 )p 0 ) ) 

[86] PlU+iNO ,NO)p (PI 41 ,P142,P143 ,P144,P145 ) 

[87] R CALCULATE THE FIRST FOUR ERROR TERMS 

[88] Elll<-Vlll*A 

[89] F[2]«-(7[2]-(BlxF[l] ))*C 

[90] F[3]-e-(7[3] ~(Pxp[i] )-(Flx£[2] ))*F 

[91] F[4]-e-(7[4] ~(CxF[l] )-(JxF[3] ))*<7 

[92] MM:N+N+1 

[93] EIN1+(RH01*ELN-1] ) + (RHOU*ElN-Ul )+7[iV] 

[94] ->(N<NOBS)/MM 

[95] E+(NOBS,l )pF 

[96] p GENERATE THE INDEPENDENT VARIABLES 

[97] XI-*- (.NOBS , 1 )p (NOBSpl ) 

[98] X2-diV0FS,l)p(W0flSWFFAWP 0 0.0625) 

[99] X+X1,X2 

[100] R GENERATE THE TRUE BETA 

[101] BET AT-*- 2 1 p 1 1 

[10 2] R GENERATE THE DEPENDENT VARIABLES 
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[103] Y*(.X+.xBETAT)+E 

[104] q LEAST SQUARES ESTIMATE OF BETA 

[105] BLS+YEX 

[106] BLSIICOUNTI1+BLSZ 1;] 

[107] BLSSZCOUNTIl+BLSl 2;] 

[10 8] n GEiVEiMZ’ff rffE i?ES JZWALS 

[109] £Mr-e-J-(X+.xSLS) 

[110] r USE FIRST ORDER P MATRIX TO TRANSFORM Y ,X , AND EH AT 

[111] r TRANSFORM Y 

[112] YSTAR1+PDW + . xj 

[113] o TRANSFORM X 

[114] XSMfllfPZW+.xX 

[115] r TRANSFORM RESIDUALS 

[116] ESTAR1+PDW+. xEHAT 

[117] o USE FOURTH ORDER P MATRIX TO TRANSFORM YSTARl , XSTARl 
[ 118] r AND ESTAR1 

[ 119] r TRANSFORM Y 

[120] YSTARn<rPn+ .xYSTARl 

[121] o TRANSFORM X 

[122] XSTAR^Pn+ .xXSTARl 

[123] r TRANSFORM RESIDUALS 

[124] ESMi?4-*-P4 + .xESMi?l 

[12 5] r CALCULATE VALUE FOR 2- STEP BETA 

[126] B2S+YSTARnmSTARU- 

[127] B2SIZCOUNTIl+B2SUil 

[128] B2SSICOUNTI1+B2SL2 ; ] 

[129] r CALCULATE THE SCHMIDT STATISTIC BEFORE 

[130] r TRANSFORMATION 

[131] DS1C+ ( (§EH AT )+. xEHAT) 

[132] DS+,(.(§EHAT) + .x(AS+.xEHAT))*DS1C 

[13 3] r COMPARE TABLED VALUE FOR SCHMIDT STATISTIC 
[134] r TO ABOVE CALCULATED VALUE FOR SCHMIDT STATISTIC 
[13 5] DEBSMT<-(DS£DSMTL ) 

[136] o IF DETECTED INCREMENT 1 

[137] DSBCNT+DSBCNT+DEBSMT 
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[138] INSBCT+ ( ( DSMTL<DS ) a ( DSMTU>DS ) ) 

[13 9] fi IF INCONCLUSIVE INCREMENT 1 

[140] STIBCT+STIBCT+INSBCT 

[141] o USE 1,4 P MATRIX TO TRANSFORM Y , X , AND EH AT 

[142] fi TRANSFORM Y 

[143] YSr4tfl4-*-P14 + .xy 

[144] fi TRANSFORM X 

[145] XSTARm+Pm + .xX 

[146] fi CALCULATE VALUE FOR 1,4 BETA 

[147] Bl4«-YS771P14gX,ST4P14 

[148] BimZCOUNTIl+Bmil;l 

[149] BmSZCOUNTI]*-Bmi2 ; ] 

[150] r CALCULATE RESIDUALS 

[151] ESTARm+YSTARlU- (XSTARin+. xfll4 ) 

[15 2] fi CALCULATE THE SCHMIDT STATISTIC AFTER 
[15 3] fi THE TRANSFORMATION 

[154] DS2D+((§ESTARin)+ .xESTARlV) 

[155] DS+, ((§ESTARm)+.x(AS+.*ESTARin))iDS2D 

[156] fi COMPARE TABLED VALUE FOR SCHMIDT STATISTIC 

[15 7] fi TO ABOVE CALCULATED VALUE FOR SCHMIDT STATISTIC 

[158] DETSMA+(DS<DSMTL) 

[159] fi IF DETECTED INCREMENT 1 

[160] DSACNT+DSACNT+DETSMA 

[161] INS ACT*- ( (DSMTL<DS )a (DSMTU>DS ) ) 

[16 2] fi IF INCONCLUSIVE INCREMENT 1 
[163] STANCT+STANCT+INSACT 

[16 4] -> (COUNTI<5 ) /INNER 

[16 5] fi CALCULATE AVERAGE VALUES FOR BETAS 

[166] BLSIAVlC0UNT01+(+/BLSI)*5 

[167] BLSSAVZC0UNT01+(+/BLSS )* 5 

[168] B2SIAVZCOUNTO]*-(+/B2SI)*5 

[169] B2SSAVZC0UNT01+(+/B2SS)*5 

[170] BmiAVZC0UNT01*-(+/BlUI)i5 

[171] BmSAVZCOUNTO] + (+/BmS)*5 

[172] +(COUNTO<VO)/ OUTER 
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[ 173 ] AF+DSACNT 

[ 174 ] BE+DSBCNT 
[17 5 ] LSI+BLSIAV 

[ 176 ] LSS<-BLSSAV 

[ 177 ] I2S<rB2SIAV 

[178] S2S<-F2SSA7 

[ 179 ] IlU+BinlAV 

[180] S14«-B14SA7 
V 
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